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Prefatory  Remarks 

This  report  is  the  first  of  a  series  of  computer  program  man¬ 
uals  arising  from  a  continuing  study  of  computer  applications  in 
the  earth  and  environmental  sciences,  including  geology,  geography, 
geophysics,  geochemistry,  and  environmental  engineering.  In  some 
of  these  subjects,  notably  oceanography,  atmospheric  science,  and 
solid-earth  geophysics,  computer  capability  has  advanced  far  beyond 
that  in  the  more  classical  fields  of  geology  and  geography,  as  well 
as  in  such  aspects  of  environmental  science  as  soil  mechanics  and 
sanitary  engineering.  Our  study  is  concerned  mainly  with  these 
more  classical  fields  in  which  computer  utilization  is  less  far 
advanced. 

Our  project  has  as  its  purposes  the  evaluation  of  developments 
in  computer  capability  in  our 'fields  through  assessment  of  present 
activities;  and  an  obligation  to  make  available  in  the  public  do¬ 
main  a  series  of  computer  programs  especially  adapted  to  the  needs 
of  workers  in  our  fields.  The  first  purpose  is  being:  met  through 
conferences  and  literature  search;  and  the  second  is  being  met  by 
reports  such  as  this,  which  will  include  programs  arising  from  our 
own  studies,  as  well  as  program,  reports  contributed  by  others 
active  in  the  field. 

The  IBM  709  program  described  here  is  based  oh  one  of  the 
earliest  programs  developed  in  the  Department  of  Geology  at  North¬ 
western  University.  It  is  almost  "classical"  in  that  it  grew  from 
an  early  machine-language  version  on  the  basic  IBM  650,  later 
modified  to  SOAP  II,  to  its  present  version  in  FORTRAN  for  the  IBM 
709.  Copies  of  the  program  in  various  stages  of  development  have 
been  widely  distributed  to  interested  workers.  Dr.  Whitten,  who 
has  developed  the  program  in  its  later  stages,  including  double¬ 
precision  computation,  has  very  kindly  prepared  this  manuscript, 
and  has  illustrated  it  with  examples  from  his  own  work  on  igneous 
rocks.  The  program  is  used  without  modification  for  studies  of 
stratigraphic,  sedimentary,  and  other  kinds  of  mappable  data.  An 
extension  of  the  program,  including  visual  trend-surface  map  output, 
is  scheduled  for  distribution  later  in  1965. 


W.  L.  Garrison 
W.  C.  Krumbein 
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A  SURFACE-PITTING  PROGRAM  SUITABLE  FOR  TESTING 
GEOLOGICAL  MODELS  WHICH  INVOLVE 
AREALLY -DISTRIBUTED  DATA 

E.  H.  Timothy  Whitten 
Geology  Department 
North-western  University 
Evanston-,  Illinois 


ABSTRACT 

The  Fortran  listing  and  program  flow  charts  are  given  in  full. 
The  method  of  utilizing  the  surface-fitting  program  is 
described  in  detail  with  the  aid  of  an  actual  example. 


INTRODUCTION 


Geological  problems  commonly  involve  large  numbers 
of  observations  which  have  been  made  at  different  map-locations. 
In  order  to  test  a  geological  response  model  on  the  basis  of 
field  observations,  an  objective  quantitative  method  must  be 
employed  to  integrate  the  data  and  to  synthesize  a  picture  of 
the  areal  variability.  Some  geological  problems  must  be  con¬ 
sidered  in  a  three-dimensional  or  poly-dimensional  framework 
(e.g.,  three  spatial  dimensions  together  with  time).  However, 
many  problems  can  be  examined  profitably  in  terms  of  a  two- 
dimensional  analysis  (or  a  series  of  two-dimensional  analyses), 
and  the  present  report  Is  restricted  to  such  cases. 

If  quantitative  measurements  (e.g.,  of  specific  gravity.) 
have  been  made  at  numerous  localities  within  a  map-area,  the 
areal  variability  can  be  expressed  visually  by  isopleths 
(i.e.,  contour-type  lines  of  equal  value).  It  is  widely 
recognized  that  such  isopleths  are  highly  subjective,  and  that 
different  draughtsmen  will  develop  wholly  dissimilar  isopleth 
maps  without  violation  of  the  observed  data.  This  situation 
could  be  rectified  by  interpolation  of  more  and  more  observa¬ 
tions;  but  in  most  practical  field  problems,  the  density  of 
sample  stations  is  prescribed  by  outside  factors  (which  are 
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Figure  1:  A  hypothetical  granite  outcrop  projected  on  to  a 
computed  mathematical  surface,  Xn  =  f(U,V).  Dots 
on  the  outcrop  area  indicate  specimen  localities 
for  which  values  of  Xn  were  observed  (After  Whitten-, 
1962,  Fig.  1). 
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usually  non-geological) .  In  consequence,  it  is  commonly  use¬ 
ful  to  obtain  an  objective  quantitative  mathematical  approxi¬ 
mation  to  the  areal  variability  shown  by  the  observations. 

The  sample  locations  can  be  defined  by  two  orthogonal 
geographic  coordinates  (independent  variables,  U  and  V). 

Ordinates  may  be  erected  .at  each  U  ,  V-location  with  heights 
proportional  to  the  quantitative  value  of  the  property 
measured  (dependent  variable,  X).  In  this  way  points  repre¬ 
senting  X  are  located  in  three-dimensions  and  the  variability 
of  X  can  be  approximated  by  a  mathematical  surface,  X=f(U,V), 
as  shown  in  Figure  1.  Various  properties  may  have  been 
measured  (e.g.,  electrical  resistivity,  quartz -content,  Pg05- 
content,  etc.),  so  that  a  family  of  surfaces  representing  X^, 

Xo,  X3,....  Xn  can  be  developed. 

Various  types  of  function  could  be  utilized  to  approxi¬ 
mate  the  best  mathematical  surface,  but  polynomials  have  been 
used  extensively  with  success.  The  present  program  uses  the 
polynomial  function: 

Xn  =  aQ  f  axU  +■  agV  f  a^  f  &4UV  +  a5v2  +  a6u3+a7Tj2v  + . (D 

To  obtain  the  coefficients  of  the  surface  which  most  closely 
approximates  the  observed  data,  the  conventional  method  of 
least  squares  is  employed.  In  practice  surfaces  of  succes¬ 
sively  higher  degree  are  computed;  the  linear  (first  degree) 


surface 

Xn  3  bQ  f  b-jU  +  b2V  ....'•  Cii) 

then  the  linear  plus  quadratic  surface  (second  degree) 

Xp,  =  Cq  -f-  c-j_U  j-  "f"  °3^  4"  c4^  +  c5^  ...... .(ill) 


and  so  on.  Successive  surfaces  of  higher  degree  (linear,  linear 
plus  quadratic,  linear  plus  quadratic  plus  cubic,....)  account 
for  larger  proportions  of  the  total  sum  of  squares  of  the 
observed  data  values  (Xn).  The  proportion  of  the  total  sum  of 
squares  accounted  for  by  a  surface  provides  some  measure  of 
its  "closeness"  or  "goodness"  of  fit. 

Commonly  it  is  convenient  to  use  the  algebraic  polynomial 
expression  to  plot  isopleths  on  the  surface  within  the  whole 
map-area.  Naturally,  the  distinction  between  such  isopleths 
and  those  drawn  with  respect  to  the  observed  data  should  be 
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kept  in  mind  clearly.  The  computed  isopleths  provide  a 
useful  method  of  assaying  the  regional  trend  inherent  in 
the  data.  In  many  geological  situations,  such  trends  are 
masked  by  the  local  variability  of  the  observed  data. 

When  the  regional  variability  has  been  approximated 
by  a  mathematical  surface,  valuable  geological  information 
can  also  frequently  be  obtained  by  noting  and  mapping  the 
departures  of  individual  observations  from  the  computed 
surface. 


DEVELOPMENT  OF  THE  METHOD  AND  TERMINOLOGY 


For  two-dimensional  graphs,  some  of  the  problems  of 
approximating  a  curve  to  a  set  of  data  points  have  long  been 
recognized.  It  is  common  practice  to  utilize  the  quantitative 
method  of  least  squares  to  develop  the  regression  line.  This 
method,  has  been  extended  to  mapped  data  and  the  technique  is 
referred  to  as  trend  surface  analysis.  Thus,  a  linear  sur¬ 
face  (Equation  il)  corresponds  to  a  regression  line  in  the 
two-dimensional  case.  The  method  of  computation  is  simple 
in  principle,  but  the  arithmetic  is  extremely  tedious  and 
was  virtually  impossible  before  the  advent,  of  high-speed  com¬ 
puting  machines. 

Oldham  and  Sutherland  (1955)  demonstrated,  the  value 
of  orthogonal  polynomials  (DeLury,  1950)  in  the  estimation 
of  regional  effects  indicated  by  mapped  data,  while  the  paper 
on  trend  surface  analysis  by  Grant ...(1957.)  may  be  considered 
definitive.  Grant  was  mainly  concerned  with  geophysical  data 
available  on  a  rectangular  grid,  but  he  showed  that  coeffi¬ 
cients  for  the  polynomial  equation  (i)  can  be  estimated  for 
irregularly-spaced  data.  The  present  program  is  designed 
for  irregularly-spaced  data  (i.e.,  observations  not  restricted 
to  rectangular  geographic  grid  intersections).  For  data 
drawn  from  sample  stations  on  a  U,V-g£id,  Grant  defined 
trend  and  residual  on  the  basis  of  a  Z'-'-array,  where  Zg 
represented  the  proportion  of  the  total  sum  of  squares  c 
of  the  mapped  variable  (Xn)  accounted  for  by  the  a^-th 
polynomial* coefficient  in  Equation  (i).  Z2 -terms  which 
contribute  significantly  to  the  total  sum  of  squares  deter¬ 
mine  the  coefficients  to  be  incorporated  in  the  complete 
trend.  The  remaining  Z2-terms  refer  to  those  coefficients 
which  comprise  the  residual.  A  typical  Z2 -array  is  shown  in 
Table  1.  It  will  be  noticed  that  the  boundary  (solid  line) 
between  the  Z^-terms  contributing  to  the  complete  trend  and 
to  the  residual  is  irregular.  Although  Grant  suggested 
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some  statistical  tests,  there  is  no  generally-accepted  pro¬ 
cedure  for  defining  this  line;  hence,  some  might  choose  a 
different  boundary  line  such,  for  example,  as  the  broken 
line  shorn  on  the  array. 

The  method  described  here  (Krumbein,  1959  )  for  com¬ 
puting  polynomial  coefficients  with  respect  to  irregularly- 
spaced  data  points  does  not  include  preparation  of  a  Z  -array, 
or  isolation  of  each  separate  coefficient.  The  linear  coeffi¬ 
cients  (Equation  ii),  the  linear  plus  quadratic  coefficients 
(Equation  iii),  the  linear  plus  quadratic  plus  cubic  coeffi¬ 
cients,  etc,,  are  computed  separately.  Since  the  complete 
trend  defined  by  Grant  (1957)  may  embrace  some,  but  not 
necessarily  all,  terms  of  several  degrees,  existing  methods 
for  irregularly-spaced  data  do  not  permit  isolation  of  the 
complete  trend. I  In  consequence,  a  surface  (e.g.,  the 
linear  plus  quadratic  surface)  is  referred  to  as  a  partial 
trend  surface  on  the  assumption  that  it  accounts  for  some  un¬ 
specified  portion  of  the  complete  trend.  When  higher  degree 
surfaces  (e.g.,  linear  through  sextic)  are  considered,  partial 

trend  surface  is  likely  to  be  a  misnomer;  the  equation  - 

probably  contains  all  the  complete  trend  terms  in  addition  to 
some  residual  terms.  Whitten  (1959)  referred  to  partial  trend 
surfaces  as  trend  components;  and  this  usage  has  been  followed 
in  numerous  publications  (e.g.,  Allen  and  Krumbein,  1962), 
although  Krumbein  (1959)  referred  to  individual  polynomial 
coefficients  as  componencs.  In  view  of  the  inadequacy  of 
partial  trend  surface,  it  was  s ingested  that  trend  component 
be  employed  in  Its  place  (Whitten,  1963). 

Whitten  (1959)  used  deviation  to  refer  to  that 
variability  not  included  in  a  partial  trend  surface  computed 
for  irregularly-spaced  data;  another  word  was  necessary  be¬ 
cause  characteristically  these  terms  are  different  from  those 
comprising  the  residual  In  Grant’s  terminology. 

Krumbein  (1959)  extended  Grant's  (1957)  method  for 
use  with  irregularly-spaced  data.  The  mathematics  are  out¬ 
lined  with  respect  to  the  computed  linear  surface 
Xn  comp.  a  bQ  -f  bpU  -f-  bgV  of  Equation  (ii). 

If  X^  is  the  observed  value  at  the  U>  ,  V n--location, 

ni  jobs  J 

then  the  deviation  at  this  point  is  (Xn.  .  -  Xn.  .  ). 

_ i3-Jobs  ijcomp _ 

^•Thls  program  has  proved  a  useful  tool  in  its  present 
form,  but  It  can  be  extended  to  determine  the  contribution 
associated  with  each  coefficient  ( cf ..  Mandelbaum,  1963). 
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To  obtain  the  least  squares  fit,  the  sum  of  squares  of  the 
deviations  must  be  minimized.  Thus, 


If 


X. 


ni,5obs  ijcompj 


\2 


biU  -  b2V 


. (iv) 


is  minimized.  It  canbe  shown  (see  Hoel,  1947,  p.  90)  that 
these  sums  of  squares  are  a  function  of  bQ,  b-j_,  and  b2  only, 
so  Equation  (iv)  can  be  expressed  as  F(bQ,  b^,  b2).  To  min¬ 
imize  P(bQ,  bp,  bg)  it  is  necessary  that 


=  3F  a  3F 

()bQ  '^b-^  2bg 

These  partial  derivatives  -are 
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W  3bg  = 

T2(X 

obs 

These  expressions  reduce 

,  b()  -  b-jU  -  b2V).(-l) 
b0  -  b-fJ  -  b2V).(-u) 

'  b0  -  V  -  b2V>-'-V> 

to  the  three  normal  equations 
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where  N  is  the  number  of  Ih  ,  V  • -locations  at  which  data  were 
obtained.  Converting  Equations  (vii)  to  matrix  form,  we  have 
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are 


Since  the 
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matrix  and  the  column,  vector 


known  from  the  original  data,.  Equation  (viii)  can  be  solved 
to  derive  the  required  coefficients  by  multiplication  of 
both  sides  by  the  inverse  of  the  [U  ,  V)  matrix.  Then: 
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The  present  FORTRAN  program  is  designed  to  compute  the  poly¬ 
nomial  linear  coefficients  by  use  of  Equation  (ix).  Similar 
equations  can  be  developed  for  coefficients  of  surfaces  of 

degree  two,  three,  four, . .  and  p.  The  present  program  is 

restricted  to  computation  of  degrees  one,  two  and  three. 

For  the  linear  plus  quadratic  plus  cubic  surface  (degree  3), 
Equation  (viii)  takes  the  form  shown  in  Equation  (x). 

It  would  appear  that  the  method  can  be  extended  to 
surfaces  of  higher  degree  provided  that  the  symmetrical  Jjj ,  V} 
matrices  can  be  inverted  satisfactorily,  or  that  other  adequate 
methods  f or  s olving  the  simultaneous  equations  are  available. 

Inversion  of  matrices  similar  to  that  in  Equation  (x) 
involve  several  problems.  In  the  course  of  extensive  testing 
the  matrix  inversion  routine  incorporated  in  this  program 
has  always  provided  adequate  inverse  matrices  when  the  compu¬ 
tation  is  effected  in  double  precision  (i.e.,  with  sixteen 
significant  figures).  The  program  can  be  extended  readily 
for  determination  of  higher  degree  polynomial  surfaces,  and 
the  existing  matrix  inversion  routine  commonly  yields  good 
results  up  to  degree  5.  For  degree  6  and  above,  an  alterna¬ 
tive  method  is  necessary  and  the  matrix-scaling  procedure 
suggested  by  Mandelbaum  (1963)  provides  one  possible  method. 

For  a  large  range  of  geological  problems,  surfaces  up  to 
degree  3  provide  adequate  trend  components.  In  fact, 
experience  has  shown  that,  in  many  cases,  use  of  surfaces  of 
higher  degree  than  three,  includes  in  the  trend  some  of  the 
variability  which  is  more  appropriately  considered  as  part  of 
the  deviation. 

To  assess  the  geological  significance  of  a  computed 
trend  component,  it  is  useful  to  develop  the  confidence 
intervals  on  the  surface.  A  method  for  determination  of  such 
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surfaces  was  described  by  Krurnbein  (1963).  por  many  pur¬ 
poses  the  sums  of  squares  associated  with  the  computed  sur¬ 
faces  provide  useful  guides  to  the  significance  of  trend  com¬ 
ponents  and  deviations.  The.  method  of  computing  these 
values  is  described  below. 


If  3<L  is  the  mean  value  of  a  set  of  observed 
nobs 

values  X„  ,  ,  then  the  deviation  from  the  mean  at  the  ii-th 
nobs 

point  is  /Xn.  .  -  X  j  ,  as  shown  in  Figure  2 A.  'The 

>  1^obs  obs/ 

sum  of  the  squares  of  all  the  deviations  is  defined  as  the 
total  sum  of  squares  of  the  observed  Xn,  i.e., 


N 

I 


(X 


n 


obs 


obs 


2 

) 


. ( xi ) 


This  expression  is  not  in  the  easiest  form  for  computation; 
it  can  be  shown  (see  Dixon  and  Massey,  1957,  p.  19)  that  Equa¬ 
tion  (xi)  is  equivalent  to: 

IT  N 

^  *nobs  -  (Z  Xnobs)  . (x 

N 


In  this  form  the  total  sum  of  squares  of  the  observed  data  can 
be  computed  very  easily. 


The  sum  of  squares 

obtained  in  a  similar  way. 

observed  values,  Xn  ,  ,  to 

obs 


of  the  computed  values,  X„  ,  is 

ncomp  - 

The  d istance  from  the  mean  of  the 
the  trend  component  (measured 


normal  to  the  U,V-plane)  at  the  ij-th  datum-point  is 
(X  -  Xn  ),  as  shown  on  Figure  23.  Hence,  the 

iJcomp  obs 

the  squares  of  the  computed  values  is: 


E(xn 

comp 


Xnobs} 


ixg  -/Zxn  \2 

comp  (  comp  J 
■  N 


sum  of 


(xiii) 


Since  X„  =  Xn  ,  the  sum  of  squares  of  the 
XIobs  ^bomp 

deviations  from  the  trend  component  (see  Figure  2C )  is  given 
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Figure  2:  Plane  with  dots  is  Xn0bs  and  plane  with-  ruled  lines  is  degree  1 
(linear)  computed  trend  component.  A  few  representative 
observed  points  (Xn)  are  shown  as  large  dots. 
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by: 


_  2 
£(Xn  -  xn  ) 
obs  corap 


£xn 

deviation 


(xiv ) 


The  proportion  of  the  total  variability  accounted  for 
by  a  trend  component  can  now  be  expressed  as  a  percentage. 

Thus,  by  using  values  derived  from  Equations  (xii)  and  (xiii), 
a  linear  surface  accounts  for: 

(Sum  of  squares  of  the  computed  Xni .  J 

- - — — - ineai —  x  100  per  cent,  (xv) 

(Total  sum  of  squares  of  the  observed  X  ) 


of  the  total  variability,  where  Xn  .  are  values  of  X  com- 

puted  for  the  linear  trend  component.  Similar  expressions  can 
be  evaluated  for  surfaces  of  each  degree  by  use  of  the  Equa¬ 
tions  (xii)  and  (xiii). 


TREND  SURFACES  IN  TESTS  OF  GEOLOGICAL  MODELS 


Within  the  physical  sciences  it  is  common  practice  to 
test  the  validity  of  conceptual  process-response  models 
with  newly-acquired  data.  This  technique  is  being  used  more 
widely  in  the  earth  sciences  (e.g..  Miller  and  Ziegler,  1958; 
Hurley,  et  al.,  1962;  Krumbein,  1962A;  Ringwood,  19 62 A,  1962B; 
Sloss,  1962;  Wyllie,  1962).  However,  few  process-response 
models  for  igneous,  sedimentary,  or  metamorphic  rocks  have 
been  evaluated  explicitly,  and  in  most  cases  the  critical 
response  and  process  factors  are  not  clearly  defined. 

In  problems  which  involve  areally  distributed  data, 
trend  surface  analysis  can  be  of  considerable  use  in  testing 
geological  models.  As  an  example  of  a  petrogenetic  model  in 
igneous  geology,  a  magmatic  granite  massif  may  be  considered. 
According  to  most  theories  a  magma  cools  inwards-  from  the 
walls  and  downwards  from  the  roof.  Such  cooling  would  tend 
to  develop  an  annular  pattern  of  mineralogical  variation  when 
the  pluton  is  examined  in  a  random  sub-horizontal  section 
(exposed  surface).  Metastable  high  temperature  mineral 
phases  would  be  preserved  where  the  magma  froze  rapidly 
(e.g.,  near  the  margins  of  the  granite).  Hence,  cafemic 
silicate  minerals  would  be  more  abundant  at  the  margins  of 
the  pluton,  and  lower  temperature  quartzo-f eldspathic 
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Figure  3:  Map  illustrating  multiple  intrusion  model 


See  text  for  explanation 
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minerals  near  the  center.  Such  fractionation  would  result 
in  both  a  chemical  zonation  and  also  a  concentric  pattern 
of  density  variation  within  the  intrusion.  Individual 
minerals  might  also  show  variation  within  the  massif.  For 
example,  more  anorthitic  plagioclase  might  be  expected 
in  the  rapidly-chilled  marginal  zones,  and  be  associated 
there  with  potassic  feldspars  with  high  t emperature  optics. 
These  phenomena  would  also  result  in  zonally-arranged 
textural  differences  within  the  pluton. 

The  process  model  for  a  cooling  magmatic  pluton  can 
be  augmented  to  include  multiple  intrusion.  Each  separately- 
intruded  unit  would  tend  to  be  annular.  The  youngest  intru¬ 
sion  would  be  expected  to  show  a  complete  sequence  of 
gradational  annular  zones,  while  the  concentric  zones  of 
older  masses  are  likely  to  be  transected  by  the  later  units. 
The  actual  contacts  between  successive  units  might  not  be 
determined  easily  in  the  field. 

Figure  3, represents  a  conceptual  model  for  the 
multiple  intrusion  hypothesis.  One  intrusion  would  result 
in  a  single  zonal  arrangement.  In  the  figure  it  is  supposed 
that  an  initial  intrusion  (1)  occurred  in  the  southeast, 
followed  by  a  younger  mass  (2)  in  the  northwest.  Both  of 
these  masses  responded  to  the  cooling  systems  by  development 
of  an  annular  pattern  of  mineralogical  zonation.  The  last 
event  represented  in  Figure  3  involved  intrusion  of  a  third 
pluton  (3),  which  transected  the  earlier  masses  (1  and  2). 
Hence,  the  last  intrusion  (3)  has  complete  zones,  whereas 
the  earlier  ones  (1  and  2)  are  incomplete. 

Other  models  could  be  erected  on  the  basis  of  dif¬ 
ferent  petrogenetic  hypotheses.  However,  as  an  illustra¬ 
tion,  the  magmatic  multiple  intrusion  conceptual  model  will 
be  examined  and  tested  with  respect  to  an  actual  example. 

The  Lacorne  granitoid  complex,  Quebec,  recently 
described  by  Dawson  and  Whitten  (1962),  affords  a  convenient 
example.  Some  recent  work  suggested  that  the  Lacorne  massif 
comprises  a.  single  intrusion,  and  some  that  it  is  a  multiple 
intrusion.  According  to  the  latter  view,  several  discrete 
units  occur  in  the  northwestern  area.  These  divergent 
views  can  be  tested  quantitatively  in  terms  of  the  con¬ 
ceptual  model  (Figure  3). 

Quantitative  modal  data  have  been  accumulated  from 
a  wide  range  of  U  ,  V-locations.  Figure  4  shows  isopleths 
drawn  to  illustrate  the  areal  variation  of  color  index 
over  the  whole  Lacorne  massif.  Commonly,  with  this  type 
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trend  c opponents  for  color  inde:c  dabs  in  the  Lacorne  "rani  bic 
nnssif,  "uehec,  Canada.  A.  linear  conponont  (accounts  for  2.; 
per  cent,  of  the  total  sun  of  squares).  D.  decree  2  coupons:!-! 
"(accounts  for  C.C  per  cent,  of  the  total  sun  of  squares).  _!• 
decree  3  component  (accounts  for  11.6  per  cent,  of  the  total 
sun  of  squares) .  (After  Cannon  and  ’..'hitton,  If 62,  'ip.  4). 
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of  data,  local  variability  tends  to  mask  the  regional 
gradients.  An  objective  method  of  screuning-out  the 
local  effects,  and  of  estimating  the  regional  trend,  is 
desirable.  If  polynomial  surfaces  are  computed  with  respect 
to  the  color  index  data,  the  trends  become  reasonably  clear. 
By  simultaneous  analysis  of  all  the  data  from  the  complex, 
the  single  intrusion  model  can  be  tested  first.  Thus, 

Figure  5  shows  the  degree  1,  degree  2,  and  degree  3  trend 
components;  these  surfaces  account  for  progressively  larger 
proportions  of  the  total  sum  of  squares  of  the  observed 
dependent  variable  (color  index).  Because  there  is  no  con¬ 
centric  pattern  of  variation.  Figure  5C  is  completely  at 
variance  with  the  model  which  involves  a  single  intrusion 
of  magma.  In  fact,  as  Dawson  and  Whitten  (1962)  pointed 
out,  this  pattern  is  anomalous  when  compared  with  the  miner- 
alogical  variation  found  in  most  other  granitic  massifs,  and 
it  suggests  multiple  intrusion. 

Recent  detailed  mapping  by  Brett  (1960)  suggested 
that  a  marked  compositional  change  occurs  in  the  area  between 
the  strongly  developed  maximum  and  minimum  illustrated  in 
Figure  5C.  Although  margins  of  separate  intrusions  were  not 
defined,  Brett's  1:12,000  map  showed  hornblende  granite 
to  the  southeast,  and  muscovite  granodiorite  to  the  northwest, 
of  the  line  of  compositional  change.  Also  a  zone  of  biotite 
granodiorite  was  mapped  along  the  junction,  but  Brett  did 
not  map  definite  contacts  or  state  whether  this  is  a  discrete 
intrusion  or  a  portion  of  either  of  the  other  two  possible 
unit  s . 


In  order  to  test  the  multiple  intrusion  model,  che 
color  index  data  from  the  entire  Lacorne  mass  were  divided 
into  two  groups  corresponding  to  the  geographic  areas 
northwest  and  southeast  of  the  zone  mapped  by  Brett,  Figure 
6A  shows  the  two  degree  3  trend  components.  Although 
eccentric,  both  areas  now  possess  a  markedly  concentric 
regional  pattern.  This  augurs  well  for  the  correctness  of 
the  model,  although  the  more  mafic  interior  of  both  areas 
is  rather  unusual  by  comparison  with  other  granites. 

Thus,  this  first  test  suggests  that  the  multiple 
intrusion  hypothesis  may  be  correct.  The  specific  geographic 
domains  which  have  been  isolated  can  be  incorporated  in  a 
slightly  more  refined  model,  which  must  be  tested  by 
quantitative  analysis  of  some  additional  variables.  Modal 
quartz  percentage  can  be  determined  easily,  and  the  trend 
components  for  this  variable  are  shown  in  Figure  6B.  These 
maps  appear  to  corroborate  the  model;  in  particular,  strong 
support  is  provided  by  the  concentric  pattern  in  the  eastern 


Figure  6:  Degree  3  trend  components  for  data  from  the 
southeastern  and  northwestern  parts  of  the 
Lacorne  massif  taken  separately.  The  percentage 
of  the  total  sum  of  squares  accounted  for  by 
the  surface  for  each  part  of  the  massif  is 
indicated  as  SS.  A.  color  index;  B.  quartz 
percentage;  C.  specific  gravity. 

(After  Dawson  and  Whitten,  1962,  Fig.  11). 
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part  of  the  Lacorne  massif.  However,  the  pattern  is  unlike 
that  of  most  other  granite  massifs  because  the  quartz  con¬ 
tent  is  lowest  in  the  center,  and  is  thus  the  opposite  of 
what  is  commonly  believed  to  be  the  case  in  most  granite 
massifs . 


The  usefulness  of  Figure  6B  as  an  independent  check 
of  the  validity  of  the  model  is  questionable  because  quartz 
percentage  and  color  index  are  both  percentage  data  which 
have  a  strong  negative  correlation  (cf.,  Chayes,  1960).  Thus, 
there  would  appear  to  be  a  strong  inherent  negative  correla¬ 
tion  between  these  two  trend  components  (Figure  6A  and  6B), 
so  that  the  quartz  percentage  surface  does  not  provide  an 
independent  reason  for  believing  that  the  model  is  correct. 

The  disharmony  between  the  trends  shown  by  (a)  these 
two  supposed  granites,  and  (b)  the  variation  which  appears 
to  be  general  in  most  other  granite  masses,  suggests  that 
additional  testing  is  necessary  before  confidence  can  be 
expressed  in  the  correctness  of  the  model.  Additional 
variables,  which  are  not  subject  to  the  closed  table 
restraints  inherent  in  the  volumetric  modal  data,  can  pro¬ 
vide  further  evidence.  Density  of  the  rocks  is  a  suitable 
variable  for  which  trend  surfaces  are  illustrated  in  Figure 
6G.  Now  there  is  no  semblance  of  a  concentric  pattern 
in  the  northwest  area,  so  it  is  doubtful  whether  a  simple 
single  intrusion  exists  there.  The  eastern  area  maintains 
a  concentric  pattern  compatible  with  that  for  the  color 
index  data  (Figure  6A)  and  does  not  compromise  the  validity 
of  the  model  for  that  area. 

With  this  information  in  hand,  careful  review  of  the 
outcrops  might  reveal  previously-unnoticed  characteristics 
which  support  the  concept  that  the  eastern  mass  is  a 
discrete  intrusion.  Because  of  the  nature  of  the  exposures, 
it  is  possible  that  the  necessary  detail  could  not  be  ob¬ 
tained.  However,  numerous  other  variables  could  be  measured 
and  subjected  to  trend  analysis.  The  variety  of  variables 
which  have  r'l ready  been  expressed  in  quantitative  terms  in 
assessing  the  areal  variability  of  different  granites 
was  reviewed  by  "/hit ten  (1963A);  other  variables  could  un¬ 
doubtedly  be  measured  to  advantage.  The  areal  variability 
of  the  temperature  of  crystallization  within  the  Lacorne 
massifs  might  provide  valuable  information  in  support  or 
contradiction  of  the  model.  Although  he  did  not  analyze 
the  results  with  trend  surfaces.  Carter  (1962)  constructed 
a  thermal  isopleth  map  for  the  Venas  granite,  Norway, 
which  indicated  the  lowest  temperatures  at  which  equilibrium 
between  co-existing  feldspars  was  attained  during  cooling  of 
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FELDS.  RATIO  DEVIATION  MAP 


GHOST  STRATIGRAPHY 
THQBR  DISTRICT,  DONEGAL 


Figure  7:  Comparison  of  palimpsestic  ghost  stratigraphy  and 
ghost  stratigraphy  in  the  ’older  granite’  of 
Donegal,  Thorr  district,  County  Donegal,  Eire. 

A.  Positive  deviations  from  the  degree  2  trend 
component  for  the  microcline  to  plagioclase  ratio 
within  granite  samples.  B.  Ghost  anticline  defined 
by  metasedimentary  rocks  enclosed  within  the 
granite  and  based  on  data  from  Pitcher  (1953A; 
1953B)  and  Pitcher  and  Read  (1959).  (After  Whitten, 
I960,  Fig.  5). 
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the  magma.  Again,  the  areal  variability  of  alkali  feldspar 
optics,  which  was  studied  in  several  Russian  granites  by 
Marfunin  (1962),  might  produce  direct  evidence  bearing  on 
the  multiple  intrusion  model  for  the  Lacorne  massif. 

The  models  discussed  above  are  susceptible  to  test 
with  trend  surface  maps  because  the  regional  variation  across 
each  unit  is  involved.  Comparable  tests  which  involve 
deviations  from  the  regional  trends  are  also  important  in 
development  of  a  complete  petrogenetic  process-response 
model;  this  is  true  for  problems  in  many  widely  different 
fields  of  geological  inquiry,  but,  initially,  the  case  of 
magmatic  granitic  intrusions  is  discussed. 

For  the  classical  model  involving  intrusion  of  a 
homogeneous  fluid.,  which  was  thoroughly  mixed  before  slow- 
cooling  and  crystallization,  the  essential  variability  would 
be  regional  (with  respect  to  the  total  intrusion).  Local  _ 
deviations  from  the  regional  areal  variability  (trend) 
would  be  essentially  random.  Deviations  might  be  expected 
to  arise  from  errors  or  inexact  measurements  during  analysis 
of  the  samples  and  from  analogous  factors  unrelated  to  the 
regional  picture.  Again,  assimilation  or  contamination  at 
the  margins  of  an  intrusion  might  result  in  local  endometa- 
morphic  modification  of  the  granite.  On  these  bases  a  petro¬ 
genetic  process-response  model  would  predict  that  the  devia¬ 
tions  are  essentially  random  over  a  major  part  of  the 
intrusion. 

Data  with  which  this  model  can  be  tested  are  sur¬ 
prisingly  sparse.  However,  Figure  7  shows  the  deviations 
from  the  d  egree  2  trend  component  for  the  ratio  of  micro- 
cline  to  plagioclase  in  the  'older  granite'  of  Donegal, 

Ireland  (Whitten,  1960).  This  map  contradicts  the  model 
because  the  deviations  show  a  distinct  pattern.  Just  prior 
to  publication  of  this  map.  Pitcher  and  Read  (1959)  established 
conclusively  that  relics  of  a  known  stratigraphic  succession 
of  Precambrian  metasedimentary  rocks  are  preserved  throughout 
this  granite  in  their  pre-granite  positions;  these  relics 
define  a  ghost  antiform.  Hence,  the  deviations  mapped  by 
Whitten  (1960)  may  be  referred  to  as  palimpsestic  ghost 
stratigraphy  and  as  a  palimpsestic  ghost  antiform,  which  are 
defined  by  local  variations  in  the  mineralogical  composition 
of  the  granitic  rock.  This  facet  of  the  quantitative 
response  model  involves  revision  of  the  proposed  conceptual 
process  model  for  this  particular  granite.  Precambrian, 
metasedimentary  rocks  which  were  where  the  granite  is  now 
must  have  radically  affected  the  composition  of  the  granite. 
Hence,  any  model  involving  intrusion  of  magma  must  involve 
tranquil  conditions  during  which  the  liquid  developed 
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inhomogeneity  prior  to  crystallization.  Alternatively, 
some  totally  different  petrogenetic  sequence  of  events 
might  have  been  involved,  in  which  case  any  magmatic  model 
would  be  proved  incorrect  by  careful  quantitative  analysis 
of  the  appropriate  variables  (response  model).  With 
financial  support  of  the  National  Science  Foundation  Grant 
G.  19633  quantitative  data  for  numerous,  additional  variables 
are  being  assembled  for  this  Irish  granite  with  a  view  to 
testing  exhaustively  various  petrogenetic  models. 

The  granite  controversy  (e.g.,  Read,  1957)  is  mainly 
involved  with  the  validity  of  rival  petrogenetic  process 
models  --  those  which  invoke  magmatic  hypotheses  and  those 
which  invoke  various  granitizational  processes.  Qualitative 
evaluation  appears  to  be  insufficient  to  resolve  the  con¬ 
troversy.  Quantitative  tests  of  specific  models  would  seem 
to  offer  a  real  opportunity  for  resolution  of  this  debate. 
Probably  for  some  granites  a  magmatic  model  will  prove  ade¬ 
quate,  but  for  others  it  may  not.  At  the  present  time  it 
is  unknown  whether  the  type  of  response  model  identified  and 
partially  described  for  the  ’older  granite*  of  Donegal  is 
unique,  or  whether  it  is  representative  of  common  relation¬ 
ships.  Deviation  maps  for  the  Malsburg  granite,  Germany, 
(Whitten,  1962)  show  that  published  K20  and  NagO  analyses 
yield  a  strongly  defined  NW.-SE.  deviation  pattern.  This  re¬ 
lationship  was  totally  unsuspected  from  the  existing  process 
models.  Hence,  current  petrological  and  petrogenetical  con¬ 
cepts  about  this  particular  granite  probably  require  revision. 
On  the  basis  of  such  a  revision,  a  new  model  could  be  con¬ 
structed,  while  the  modal  and  chemical  data  recently  published 
by  Leible  (1959),  Mehnert  (I960),  Mehnert  and  Willgallis 
(1961),  and  Rein  (1961)  would  form  an  initial  basis  for  quan¬ 
titative  tests  of  the  new  model. 

Wadsworth  (1963)  studied  the  areal  variability  of 
textural  variables  in  the  Twelve-foot  Fall  quartz  diorite, 
Wisconsin,  with  the  aid  of  trend  surface  analysis.  Work  with 
many  additional  types  of  variables  must  supplement  analyses 
of  the  standard  modal  and  chemical  measurements  if  an  ade¬ 
quate  response  model  is  to  be  constructed. 

At  present  many  qualitative  and  subjective  judgments 
about  the  composition  and  variability  of  rock  masses  have  to 
be  relied  upon  in  petrogenetic  theory.  If  adequate  data  can 
be  made  available,  trend  surface  and  deviation  maps  provide 
objective  estimates  which  are  valuable  as  firm  foundations 
for  evaluation  of  petrogenetic  models  for  granite  bodies. 

Although,  this  discussion  has  been  illustrated  by 
granites,  it  must  be  emphasized  that  trend  surface  and 
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deviation  maps  have  value  in  many  other  fields  of  the  earth 
sciences.  The  earliest  deviation  maps  were  published  for 
stratigraphic  facies  maps  by  Krumbein  (1956,  1968B),  Grant 
(1957)  showed  that  positive  deviations  based  on  observed 
Bouguer  anomalies  are  associated  with  sub-surface  sulphide 
mineralization.  On  the  basis  of  deviation  maps  prepared 
for  the  heavy  mineral  content  of  the  top  Ashdown  Pebble  Bed, 
England,  Allen  and  Krumbein  (1962)  developed  a  paleogeo- 
graphic  reconstruction  (process  model)  for  the  Wealden  area. 

Most  geological  problems  are  susceptible  to  quantita¬ 
tive  analysis.  Commonly,  the  philosophy  of  objectives  in 
geological  mapping  is  not  clear.  As  Chayes  and  Suzuki  (1963) 
suggested,  greater  clarity  could  result  from  closer  definition 
of  the  objective  ( cf Whitten,  1963B).  In  many  cases  the 
objective  will  remain  obscure  until  the  processes  involved, 
and  the  responses  to  these  processes,  are  expressed  quanti¬ 
tatively  in  the  context  of  a  unifying  model.  The  model  will 
have  to  be  modified  as  more  is  learned  of  the  processes  and 
the  responses  associated  with  particular  geological  problems. 
In  studies  of  present-day  sedimentation,  both  the  processes 
and  the  responses  can  commonly  be  analysed.  The  processes 
involved  in  the  formation  of  granites  are  essentially  con¬ 
jectural.  In  the  erection  of  process-response  models  for  the 
genesis  of  granites  magmatic,  granitizational,  or  some  other 
processes  may  be  involved.  The  validity  of  such  models  has 
to  be  assessed  solely  on  the  results  of  quantitative  analysis 
of  the  many  attributes  of  the  response  products  (the  granitic 
rocks ) . 


The  present  program  provides  a  useful  method  of  eval¬ 
uating  the  quantitative  behavior  of  response  characteristics, 
and  hence,  of  testing  petrogenetic  and  other  geologic  models. 


DEVELOPMENT  OF  THE  FORTRAN  PROGRAMS 


Professors  W.  C.  Krumbein  and  D.  Harris  prepared  a 
surface-fitting  program  in  MACHINE  LANGUAGE  for  the  Basic 
IBM  650  in  1957.  This  program  would  compute  degree  1,  2,  and 
3  trend  components  with  respect  to  four  dependent  variables, 
and  it  effected  the  calculations  in  single  precision.  This 
program  was  rewritten  in  SOAP  and  filed  in  the  IBM  London 
Library  by  Krumbein  and  C.  Faulkner  in  1960  under  file  number' 
60705,  and  was  subsequently  catalogued  under  the  number 
8.3.001  in  IBM  New  Y0rk  Library. 
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In  June,  1961,  Mr.  R.  Axelrod  began  re-writing  the 
program  for  Krumbein  in  FORI  RAN  language  for  use  on  the 
IBM  709.  By  February,  1962,  Whitten  had  converted  the 
entire  FORTRAN  program  to  double  precision  and  extended 
it  for  simultaneous  use  with  eight  variables.  In  a ddi- 
tion,  Whitten  prepared  the  program  designated  as  Part  II 
in  this  report  for  use  in  conjunction  with  the  main  surface¬ 
fitting  program. 

Apart  from  these  main  steps  in  the  evolution  of  this 
program,  minor  revisions  have  been  made  continuously.  Through¬ 
out  this  work  the  staff  of  the  Northwestern  University  Com¬ 
puting  Center  have  given  much  advice;  special  mention  should 
be  made  of  assistance  given  by  Mrs.  0.0.  Benson  and  Mr.  H. 

Rymer  of  Northwestern  University,  and  by  Miss  S.  Hitz  and 
Dr.  C.  Faulkner  of  IBM. 

Financial  support  for  this  work  has  been  received 
from  the  Graduate  School  of  Northwestern  University,  the 
National  Science  Foundation  (grant  to  Whitten,  Project  Num¬ 
ber  G-19633),  and  the  Office  of  Naval  Research  (Contract 
Nonr-1228(26) ;  Task  no.  389-135).  This  continuing  assistance 
is  gratefully  acknowledged. 


OUTLINE  OF  THE  PROGRAM  AND  CAPACITY 


The  present  program  is  designed  to  operate  with  one 
to  eight  sets  of  mapped  variables  (dependent  variables). 
There  is  no  upper  limit  to  the  number  of  data-points 
which  can  be  utilized.  The  lower  limit  to  the  number  of 
data-points  is  prescribed  by  the  distribution  of  the  data- 
points  and  the  precision  with  which  the  surface  must  be 
defined,.  However,  as  an  empirical  guide,  at  least  three 
times  the  number  of  data-points  as  the  number  of  coeffi¬ 
cients  should  be  used. 

Computation  of  the  polynomial  coefficients  is 
facilitated  if  the  axes  are  oriented  so  as  to  reduce  un¬ 
intentional  alignment  of  the  U  and  V  values  of  the  data- 
points.  Commonly,  the  origin  Is  placed  at  the  top  left- 
hand  corner  of  the  map,  so  that  U  increases  downwards  and 
V  from  left  to  right.  This  choice  was  based  on  preserva¬ 
tion  of  a  sense  of  direction  in  the  map  similar  to  that 
U3ed  for  g ridded  data,  which  are  commonly  analysed  as  rows 
and  columns  starting  from  the  upper  left  (Krumbein,  1960, 
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p.  360),  It  is  advantageous  to  maintain  this  orientation., 
although  others  can  be  found  in  the  literature. 

The  FORTRAN  program  can  be  considered  conveniently 
as  a  series  of  sections,  thus: 


PART  I 

Phase  1 

Phase  2 

Phase  3 

Phase  4 


Phase  5 


Preparation  of  the  |jj  ,  v]  matrix  and 
the  column  vector  [xn]  . 

Inversion  of  the  degree  1,  2,  and  3 
['ll  ,  vj  matrices,, 

Multiplication  of  the  inverse  matrices 
and  column  vectors  [xn]  to  obtain  the 
degree  1,  2,  and  3  coefficients. 

Use  of  the  coefficients  to  compute 
values  of  Xn  at  each  datum-point  and 
to  calculate  the  deviations  at  each  of 
these  map -points. 

Preparation  of  the  sums  of  squares  sum¬ 
mary  for  the  surfaces  of  each  degree. 


PART  II 

Use  of  the  coefficients  to  calculate  a 
network  of  computed  values  over  the  entire 
map -are a. 


The  actual  FORTRAN  program  is  listed  in  Appendix  I,  and 
the  flow  charts  are  shown  in  Figures  8  and  9, 
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It  will  be  noticed  that  most  of  the  operations  in 
Part  I  are  completed  in  double  precision,  which  insures  use 
of  sixteen  significant  figures.  This  is  necessary  because 
use  of  eight  significant  figures  only  (single  precision)  in¬ 
troduces  serious  rounding  errors  at  several  stages  of  the 
program. 


In  practice  Part  II  is  an  invaluable  adjunct  to  the 
main  surface-fitting  program  (Part  I)  because  commonly  it  is 
difficult  to  draw  isopleths  accurately  on  the  basis  of  values 
computed  at  the  observed  data-locations  only.  Isopleths  for 
the  computed  surfaces  are  mathematical  curves.  These  are 
drafted  more  easily  if  a  network  of  computed  values  is  laid 
over  an  area  somewhat  greater  than  that  from  which  the  ob¬ 
served  data  were  collected.  With  Part  II,  values  can  be  com¬ 
puted  for  a  rectangular  grid  of  up  to  35  points  in  both  the 
U  and  the  V  directions.  The  simple  calculations  involved  in 
Part  II  do  not  require  double  precision. 


PREPARATION  TO  RON  THE  PROGRAM 
(PART  I) 


Data  cards: 

A  separate  data  card  is  used  for  each  sample  location; 
the  U,V -co ordinates  and  the  data  for  one  through  eight  dependent 
variables  recorded  at  each  locality  are  included  on  the  same 
80-column  IBM  card.  The  format  requires: 


Columns 


1 

through 

6 

7 

?i 

10 

11 

!! 

16 

17 

tr 

22 

23 

?! 

28 

29 

34 

-  Project  number  (either  numeric  or 
alphabetic) 

-  Sample  point  identification  (either 
numeric  or  alphabetic) 

-  Geographic  coordinate  U  (an  independent 
variable ) 

-  Geographic  coordinate  V  (an  independent 
variable ) 

-  Dependent  variable  Xq 

-  Dependent  variable  Xg 


Succeeding  groups  of  six  columns  accommodate  X3  through 
Xg.  Columns  71  through  80  are  left  blank. 

A  typical  card  with  eight  dependent  variables  (all 
identical  for  illustration)  w?uld  have  the  following  appearance: 
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.  5000 1 , 5000 1 , 5000 1 , 5000 1 ,  <00  i . !;' 

I  I  I  I  I  I 
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1  l  3  4  5  S  7  e  S  10  II  12  13  14  15  It  17  II  IS  20  71  2!  23  24  25  2$  27  21  23  30  31  32  33  34  35  35  37  31  33  40  41  42  43  44  45  46  47  41  49  SO  51  52  53  54  55  56  57  58  59  60  01  62  51  64  65  65  57  65  69  70  71  72  73  74  75  76  77  78  79  50 
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2  2  2  2  2  2  2  2  2|  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2 

|  3333 33333 33 3|33333|33333|33333|333 3 3|333 3 3|3 3 333|3 3 3 33(33333 |33333 | 333 33333 33333  | 

4(  444444  ♦  444444444  4' 4-4  44  4. 4444444444444444444  +  44444  4' 4444444444  4' 4444444444444444444  S 
5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5|  5  5  5  5  5|  5  5  5  5  5|  5  5  5  5  5|  5  5  5  5  5|  5  5  5  5  5|  5  5  5  5  5J  5  5  5  5  5|  5  5  5  5  5  5  5  5  5  5  5  5 
06666666666666666680668680666666066666666666666666660666666666666666666666666666 
77777|77777777777777777773777777777777/77777777777777777777777777777777777777777 
|  8  8  8  8  S  8  8  8  8  8  8  8|  8  8  8  8  8|  8  8  8  8  8|  8  8  8  8  8|  8  8  8  8  S|  8  8  8  8  8|  8  8  8  8  6|  8  8  8  8  8|  8  8  8  8  8|  8  8  8  8  8  |  8  8  8  8  8  8  8  8  8  8  8  8  8 

99999999 9999 9 99 9999999999999999999999909 999999999 9099903 99999999 99999999999 99 999 

J  2  3  4  5  C  7  8  9  10  11  12  13  14  IS  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  V  ^75  76  77  78  79  80 


Each  variable  should  be  expressed  to  the  same  number 
of  decimal  places.  It  is  often  convenient  to  punch  the 
decimal  on  the  data  cards  but  this  is  optional  as  the  posi¬ 
tion  of  the  decimal  is  specified  on  the  second  master  card. 
The  program  cannot  accommodate  'no  data,'  because  a  blank 
space  is  interpreted  as  a  zero  observation.  The  program 
sets  no  limit  to  the  number  of  data  cards  (and  thus  to  the 
number  of  U,V-locations ) . 

If  less  than  eight  dependent  variables  are  involved, 
more  space  can  be  used  on  the  data  card  for  each  variable; 
although  the  format  of  columns  1  through  10,  and  71  through 
80  remain  unchanged.  For  example,  if  four  dependent  variables 
are  involved,  each,  variable  could  be  assigned  to  10-digit 
words,  so  that  U  occupies  columns  11-20,  V  columns  21-30, 
etc . 


The  dab  a  de ck ; 


The  data  deck  is  assembled  in  the  following  manner j 
1 .  *DATA 


2 


Four  TITLE  CARDS 
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3.  Two  MASTER  CARDS 

4.  Deck  of  data  cards 

5.  One  WIRES  CARD 

6.  Deck  of  data  cards  (exact  duplicate  set  of  4  above) 

7.  One  NINES  CARD  (exact  duplicate  of  5  above) 

If  more  than  one  problem  is  to  be  executed  with  this 
program  in  the  same  run,  successive  decks  may  be  assembled 
according  to  2  through  7  above,  and  be  placed  behind  7  of 
the  first  deck.  Any  number  of  successive  problems  can  be 
executed  in  this  manner. 

It  will  be  noticed  that  two  identical  data  decks 
(4  to  6  above)  are  required..  This  permits  an  unlimited  num¬ 
ber  of  geographic  locations  to  be  employed,  and  this  benefit 
outweighs  the  slight  inconvenience  of  duplicating  the  data 
cards.  The  program  could  be  modified  to  store  the  data  card 
information  ready  for  second  reading  (B  in  Plow  Sheet,  Figure 
1);  this  would  obviate  the  need  for  a  second  deck  of  data 
cards,  but  involve  a  limit  on  the  number  of  data  cards  which 
could  be  accommodated. 


Title  cards: 

The  format  is  (9  A6)  so  numeric  and  alphabetic 
characters  can  be  utilized;  typical  cards  might  have  the 
following  form; 

1  WHITTEN  PROJECT  040060 

0  MODAL  DATA  FOR  LACORNE  PLUTON,  QUEBEC 

0  XI  =  QTZ;  X2  =  COLOR  INDEX;  X3  -  FELD.  RATIO 

0  LQC -TREND  SURFACES,  N  =  48 


Master  cards; 

The  first  card  has  the  format  (A6,  A4,  II,  FI. 0)' and 
must  provide  the  following  information; 
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Columns  1  through  10  -  project  control  number  (numeric  or 

alphabetic  characters)  identified  as 
NPR0J1  and  NPRQJ2 

Column  11  -  Number  of  dependent  variables  (1 

through  8)  identified  as  NX 

Column  12  -  Insert  1  if  require  coefficient  cards 

to  be  punched  by  program  in  addition 
to  list  of  coefficients  printed  auto¬ 
matically.  If  cards  are  not  required, 
leave  column  blank.  Identified  in 
program  as  PC ECO. 

The  second  master  card  has  a  format  of  (10X,  10A6), 
fmd  the  statement  is  identified  in  the  program  as  VPT  (variable 
format  statement).  The  card  is  prepared  as  follows: 

Columns  1  through  10'  -  Project  control  number  (numeric  or 

alphabetic  characters)  not  used  in  pro¬ 
gram,  but  used  to  identify  this  card. 

Columns  11  on  -  Statement  of  the  format  employed  for 

the  data  cards.  A  typical  statement 
would  be  ( 6XA4, 10F6. 3), which  implies 
that  (i)  columns  1  through  6  of  data 
card  are  not  read,  (ii)  columns  7 
through  10  are  read  and  contain  numeric 
or  alphabetic  sample  point  identifica¬ 
tion,  and  (iii)  columns  11  onwards  con¬ 
tain  10  words  of  6  letters  all  expressed 
to  3  decimal  places  (this  would  be  ap¬ 
propriate  for  U,  V,  and  eight  Xn).  With 
U,  V ,  and  two  Xn  entered  to  two  decimal 
places  in  10-letter  words,  the  state¬ 
ment  would  be  ( 6XA4, 4F10.2 ) . 

Thus.,  a  typical  set  of  two  master  cards  might  be 
040060000041  (i,e.,  four  variables  and  coefficient 

0400600000( 6XA4, 10P6.3 )  cards  wanted) 

( N o t e :  the  format  statement  specifies  the  same  number  of 
decimal  places  for  the  variables  U,  V  and  Xn). 


Nines  card: 

The  program,  searches  for  the  value  99999.  to  indicate 
that  all  data  cards  have  been  read.  When  U,  V,  and  Xn  are 
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expressed  as  six-letter  words,  the  nines  card  requires  99999. 
in  columns  11  through  16.  If  ten-letter  words  are  employed 
000099999.  is  punched  in  columns  11  through  20,  and  similarly 
for  any  other  word  length  99999.  is  punched  in  the  U-position 
of  the  data  cards. 

In  this  FORTRAN  program  all  input  is  read  from  Tape 
5.  Output  is  put  on  Tape  6  when  printed  listing  is  required 
and  on  Tape  7  when  card  output  is  requested. 


OPERATION  OF  THE  PROGRAM 


(PART  I) 


The  sequence  of  operations  can  be  read  from  the  flow 
sheet  (Figure  8).  This  sequence  and  the  output  of  the  program 
is  outlined  below. 


For  purposes  of  illustration,  identical  values  are 
used  for  X.±,  X2.»-».Xq  at  each  datum-location  in  the  following 
example.  This  has  been  done  solely  because  the  illustrative 
output  is  more  easy  to  read,  but  in  practice  the  several  X_ 
may  be  completely  independent,  although  each  is  a  dependent 
variable  with  respect  to  U  and  V. 


The  four  title  cards  are  read  first  and  this  informa¬ 
tion  is  immediately  printed  out  as  shown  in  Table  2. 


Following  some  initialization,  the  information  con¬ 
tained  on  the  two  master  cards  is  read  and  stored  for  use 
throughout  the  program  (point  A  on  flow  sheet.  Figure  8). 

The  program  now  begins  to  read  the  first  deck  of 
data  cards  and  assembles  the  [U,V.]  matrix  and  the  column 
vector  [Xn]  .  During  initialization  NPHAS  was  set  equal  to 
1,  so  that  as  each  data  card  is  read  the  information  is 
immediately  printed  out  on  one  line.  The  powers  and  products 
of  U  and  V  for  this  datum -location  are  then  computed  (37 
in  program  listing)  and  stored  in  the  10  x  10  matrix  P(I,J) 

(point  C  of  flow  sheet.  Figure  1)  and  the  column  vector 
3(1, J).  This  column  vector  [X  ]is  actually  a  matrix  of  size 
10  x  NX,  where  NX  is  defined  by  the  number  of  dependent 
variables,  so  the  matrix  actually  stores  one  column  vector- 
for  each  dependent  variable  (Xn). 

Then  the  second  data  card  is  read.;  the  additional 
information  is  added  to  that  already  in  the  ?(I,J)and  B(  I,  J«)matriees 
Successive  data  cards  are  dealt  with  until  the  nines  card  is 
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TABLE  2 


WHITTEN  PROJECT  0401070000 

PROBLEM  TO  ILLUSTRATE  PROGRAM 
XI  THRU  X8  ALL  IDENTICAL  SYNTHETIC 

LQC  DOUBLE  PRECISION  8-  VARS.  ‘ 


TABLE  3 


DATA 


CONTROL 

u 

0001 

0. 

0002 

1.000 

0003 

2.000 

0004 

3.000 

0005 

0. 

0006 

1.000 

0007 

2.000 

0008 

3.000 

0009 

0. 

0010 

1 . 000 

0011 

2.000 

0012 

3.000 

0013 

0. 

0014 

1.000 

0015 

2.000 

0016 

3.000 

0017 

0. 

0018 

1.000 

0019 

2.000 

0020 

3.000 

V 

XI 

0. 

1.000 

0. 

1.500 

0. 

1.000 

0. 

2.000 

1.000 

2.000 

1.000 

2.500 

1.000 

3.000 

1.000 

2.000 

2.000 

2.000 

2.000 

2.500 

2.000 

3.000 

2.000 

2.000 

3.000 

0.020 

3.000 

0.500 

3.000 

2.000 

3.000 

1.000 

4.000 

0. 

4.000 

0.500 

4.000 

1.000 

4.000 

0.500 

X2 

X3 

1.000 

1.000 

1.500 

1.500 

1.000 

1.000 

2.000 

2.000 

2.000 

2.000 

2.500 

2.500 

3.000 

3.000 

2.000 

2.000 

2.000 

2.000 

2.500 

2.500 

3.000 

3.000 

2.000 

2.000 

0.020 

0.020 

0.500 

0.500 

2.000 

2.000 

1.000 

1.000 

0. 

0. 

0.500 

0.500 

1.000 

1.000 

0.500 

0.500 

X4 

X  5 

1.000 

1.000 

1.500 

1.500 

1.000 

1.000 

2.000 

2.000 

2.000 

2.000 

2.500 

2.500 

3.000 

3.000 

2.000 

2.000 

2.000 

2.000 

2.500 

2.500 

3.000 

3.000 

2.000 

2.000 

0.020 

0.020 

0.500 

0.500 

2.000 

2.000 

1.000 

1.000 

0. 

0. 

0.500 

0.500 

1.000 

1.000 

0.500 

0.500 

X6 

X7 

1.000 

1.000 

1.500 

1.500 

1.000 

1.000 

2.000 

2.000 

2.000 

2.000 

2.500 

2.500 

3.000 

3.000 

2.000 

2.000 

2.000 

2.000 

2.500 

2.500 

3.000 

3.000 

2.000 

2.000 

0.020 

0.020 

0.500 

0.500 

2.000 

2.000 

1.000 

1.000 

0. 

0. 

0.500 

0.500 

1.000 

1.000 

0.500 

0.500 

TABLE  4 


UV  CUBIC 

l 

0. 20000000E 
0.  1 2000000E 

02 

03 

0. 30000000E 

0.  18000000b 

02 

03 

0 . 400000006 
0.  139999996 

02 

03 

0.700000006 

0.180000006 

02 

03 

0. 59999999E 
0,400000006 

02 

03 

2 

0. 30000000E 
0. 18000000E 

02 

03 

0. 7  OOOOOOOE 
0.48999999b 

02 

03 

0.599999996 

0.360000006 

02 

03 

0.180000006 

0.419999996 

03 

03 

0,  13999999E 
0.60000000b 

03 

03 

3 

0. 40000000b 
0. 40000000b 

02 

03 

0. 59999999E 
0. 36000000E 

02 

03 

0.  12000000E 
0.619999996 

03 

03 

0.139999996 

0.600000006 

03 

03 

0.  1 8000000E 
0.  14 160000E 

03 

04 

4 

0.70000000b 
0.4 1999999b 

02 

03 

0. 18000000E 
0. 13800000b 

03 

04 

0. 139999996 
0. 98000000E 

03 

03 

0. 48999999E 
0. 108000006 

03 

04 

0.36000000E 
0.  1  3999999E 

03 

04 

5 

0. 59999999E 
0. 60000000E 

02 

03 

0.13999999b 
0. 9R000000L 

03 

03 

0. 130000006 
0.108000006 

03 

04 

0.360000006 
0..  139999996 

03 

04 

0.4199 9 99 9G 
0.2  1240000E 

03 

04 

6 

0. 12000000E 
0. 14 I60000L 

03 

04 

0. 18000000E 
0. 1 0800000E 

03 

04 

0.400000006 
0. 139999996 

03 

04 

0.419999996 

0.212400006 

03 

04 

0.60000000E 

0.52000000b 

03 

04 

7 

0. 18000000b 
0. 10800000C 

03 

04 

0 • 48999999E 
0. 39700000C 

03 

04 

0.360000006 

0.276000006 

03 

04 

0.138000006 

0.29400000E 

04 

04 

O'. 98000000E 
0. 36000000E 

03 

04 

8 

0.  13999999b 
U.  13999999C 

03 

04 

0. 36000000E 
0.27600000b 

03 

04 

0.419999996 

0.294000006 

03 

04 

0.980000006 

0.360000006 

03 

04 

0. 10800000b 
0.4'9560000t 

04 

04 

9 

0.  18000000b 
0.2  1 2 400 00 L 

03 

04 

0.41999999b 

0.29400000C 

03 

04 

0.600000006 

0 • 36000000  E 

03 

04 

0.  10800000E 
0.495600006 

04 

04 

0. L3999999E 
0.78000000b 

04 

04 

10 

0. 40000000L 
C.52000000F 

03 

04 

0. 60000000C 
0. 36000000G 

03 

04 

0. 141600006 

0 .493600006 

04 

04 

0.139999996 

0.780000006 

04 

04 

0.212400001- 

0.19560000b 

04 

05 

X8 

1,000 
i  •  500 
I. 000 
2.000 
2.000 
2.500 
3.000 
2.000 
2.000 
2.500 
3.000 
2.000 
0.020 
0.500 
2.000 
1.000 
0. 

0.500 
I. 000 
0.500 
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TABLE  5 


VECTOR  X(S) 


L 

0*300199996 

02 

0.30019999b 

02 

0.300199996 

02 

0.300199996 

02 

0.30019999b 

02 

0.30019999E 

02 

0.300199996 

02 

0.30.0199996 

02 

2 

0.50000000E 

02 

0. 50000000E 

02 

0. 50000000E 

02 

0. 50000000E 

02 

0.500000006 

02 

0.500000006 

02 

0.500000006 

02 

0 • 50000000E 

02 

3 

0.470599996 

02 

0. 4  7059999E 

02 

0.47059999E 

02 

0.470599996 

02 

0. 47059999E 

02 

0.47059999E 

02 

0.470599996 

02 

0.470599996 

02 

4 

0. 114999996 

03 

0.1 149Q999E 

03 

0.1 1499999E 

03 

0.114999996 

03 

0. 1.14999996 

03 

0.  U499999E 

03 

0. 1 14999996 

03 

0. 11499999E 

03 

5 

0.820000006 

02 

0.820000006 

02 

0.820000006 

02 

0.820000006 

02 

0.820000006 

02 

0.820000006 

02 

0.820000006 

02 

0.82000000E 

02 

6 

0. Ill 17999E 

03 

0.111179996 

03 

0.111179996 

03 

0.11  117999E 

03 

0. 11  H79996 

03 

0.111179996 

03 

0. 111179996 

03 

0. 111179996 

03 

7 

0.290000006 

03 

0.290000006 

03 

0.290000006 

03 

0.290000006 

03 

0 . 29000000E 

03 

0.290000006 

03 

0. 29000000E 

03 

0. 29000000E 

03 

8 

0. 186000006 

03 

0. 18600000E 

03 

0.186000006 

03 

0.186000006 

03 

0.186000006 

03 

0. 186000006 

03 

0.186000006 

03 

0. 18600000E 

03 

9 

0.204000006 

03 

0. 20400000E 

03 

0.204000006 

03 

0.204000006 

03 

0.204.000006 

03 

0.204000006 

03 

0.204000006 

03 

0 . 20400000E 

03 

10 

0.308539996 

03 

0.308539996 

03 

0.30853999E 

03 

0.308539996 

03 

0.308539996 

03 

0.308539996 

03 

0  .'3  08  53  99  9  6 

03 

0 • 30853999E 

03 

TABLE  6 

INVERTED  MATRIX 

OF  DEGREE  l 

1 

0.239999996-00 

-0.599999996-01 

-0.49999999t-0l 

2 

-0.599999996-01 

0.399999996-01 

-0. 

3 

-0.499999996-01 

-0. 

0.24999999E-01 

TABLE  7 

INVERTED  pihIkia  OF 

DEGREE  2 

1 

0.54L428576  00 
0.357142856-01 

-0.330000006-00 

-0 • 202857 14E-00 

0.499999996-01 

0.59999999C-01 

2 

-0.330000006-00 

-0. 

0. 56999999E  00 

0. 59999999E-01 

-0.150000006-00 

-0.399999996-01 

3 

-0. 282857 L4E-00 
-0.7L4285716-01 

0.599999996-01 

0.35571428E-00 

-0. 

-0.300000006-01 

4 

0.499999996-01 

0. 

-0.15000000C-00 

-0. 

0.499999996-01 

0. 

5 

0.599999996-01. 

0. 

-0.399999996-01 

-0.300000006-01 

0. 

0.200000006-01 

6 

0. 357142856-01 

0. 178571436-01 

-0. 

-0.7142857U-01 

0. 

0. 
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reached  at  the  end  of  the  data  deck,  by  which  time  a  complete 
list  of  the  data  cards  has.  been  made  in  the  format  shown  in 
Table  3. 


Immediately  the  nines  card  is  read,  the  program  calls 
for  a  print  out  of  P(I,J),  which  comprises  the  F_U,V]  degree 
3  matrix  (see  Equation  x).  This  information  is  listed 

in  floating  point  with  the  first  row  of  the  matrix  across 
the  page  in  two  lines  identified  by  1.  Subsequent  rows  are 
identified  by  2,  3,  4, .....10  (see  Table  4). 

The  column  vectors  for  each  Xn  are  then  printed  in 
floating  point  (Table  5).  The  ten  entries  of  each  column 
are  listed  down  the  page  (identified  by  1,  2,  3, .....10), 
and  successive  columns  across  the  page  correspond  to  Xp,  X2, 

X3, . X8  (those  for  Xg  through  Xq  being  listed  on  a  second 

line ) . 


Initialization  for  Phase  2  of  Part  I  is  then  com¬ 
pleted  and  the  program  reads  control  statements  in  preparation 
for  inversion  of  the  linear  3x3  Cu,Vj  matrix  (D  in  flow 
sheet).  The  linear  jjJ, V]  matrix  comprises  the  top  left 

3x3  block  of  the  10  x  10  matrix  stored  in  P(I,J).  The  in¬ 
verted  3x3  matrix,  called  A(I,J)  in  the  program,  is  printed 
out  in  floating  point  as  shown  in.  Table  6.  This  inverted 
matrix,  A(I,J),  is  stored  in  Y( I, J,NDEG) ,  in  which  NDEG  is 
defined  as  1  (degree  1). 

The  program  proceeds  to  802  and  reads  the  degree  2 
controls  before  returning  to  the  matrix  inversion  routine 
(D  on  the  flow  chart)  to  invert  the  degree  2  matrix,  i.e., 
the  top  left  6x6  section  of  the  main  matrix, 

P(I,J).  The  Inverted  6x6  matrix,  A(I, J),  is  printed  out 
in  floating  point  (see  Table  7)  and  stored  in  Y(I,J,NDEG) 
with  the  degree,  NDEG,  defined  as  2.  The  sixth  column  of 
A(I,J)  is  listed  on  a.  lower  line. 

Proceeding  to  803,  the  degree  3  controls  (including 
NDEG  -  3)  are  read  and  the  whole  10  x  10  P(I,j)  matrix  is 
inverted.  The  inverse  A(I,J)  is  printed  out  and  stored  in 
Y(  I,  J,NDEG) .  Each  row  of  A(.I,J)  requires  two  lines  to 
print  out  as  shown  in  Table  8.  Because  P(I,J)  is  always 
symmetrical  about  the  principal  diagonal,,  the  degree  1,  2, 
and  3  inverse  matrices,  A(I,J),  should  also  be  symmetrical 
about  their  principal  diagonal.  The  lists  of  the  three 
inverted  matrices  enable  the  symmetry  of  each,  matrix  to  be 
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1 

2 

3 

4 

5 

6 

7 

8 

9 


10 


1 

2 


3 


1 

2 

3 

4 

5 

6 


TABLE  8 


INVER  I  EH  M A  I  K I  X  I'JF  DEGREE  3 


U. 80500000b  00 

0. 22499999b-00 

-0.872380951  00 
-0.33333333C-01 

-0.769166668  00 
-0.499999998-01 

0.300000008-00 
-0. 428571438-01 

-0.87238095b  00 
-0.428571431-01 

0.398158721:  01 
0.522222228  00 

0 . 38 142857E-00 

0. 130000008-00 

-0.279999991  01 
0.285714281-01 

-0.769166661  00 
-0.1095833)1  01 

0. 381428578-00 
0.292112011-15 

0.217902771  01 
0.249999998-01 

-0.499999991-01 

0.857142851-01 

U. 30000000L-00 
-0.296039478-18 

-0.27999999E  01 
-0.500000008  00 

-0.499999998-01 

-0.499999998-01 

0.239999991  01 

0. 

0.381423578-00 

0.857142858-01 

-0.604285711  00 
-0.655278308-15 

-0.44785714E-00 
-0.749999998-0 1 

0.  150000008-00 
-0.571428576-01 

0. 224999991-00 

0. 6 7500000E  00 

-0.428571438-01 

0.300518438-31 

-C .  10958  3338  01 

0. 

-0.135233291-30 

-0.214285718-01 

-0. 33333333E-01 

0. 

0.522222228  00 

o.iiiniioi-oo 

-0. 

0. 

-0.499999991-00 

-0. 

-0.499999998-01 

0.  1.48029731-  1 5 

0. 15000000E-00 

0. 1 97372988-15 

0.249999991-01 

0.249999991-01 

-0.499999991-01 

-0. 

-0.42857143L-01 

-0.214285718-01 

0.285714281-01 

-0.200345621-31 

0.857142851-01 

-0. 

0.901555318-31 
0. 142H5714E-01 

-0.208333338-01 
-0.  1041667, 6E-00 

0. 

-o. 

0.  149305558-00 

-0. 

0. 

0. 

TABLE  9 

CHIEFS  DEGREE  1 

0.185179998  01 
0.185179998  01 

0. 185179996  01 
0.185179998  01 

0.185179991  01 
0.183179996  01 

0.183179998  01 

0. 198800008-00 

0. 198800008-00 

0. 198800001-00 

0. 198800008-00 

0. 19B80000E-00 

0.  198000001-00 

0. L9880000E-00 

-0.324499991-00 

-0.324499998-00 

-0. 32449999E-00 
-0.324499998-00 

-0.324499991-00 

-0.324499991-00 

-0.324499991-00 

. 

TABLE  10 

CUEFFS  DEGREE  2 

0.108314288  01 
0.108314281  01 

0.108314288  01 

0. 108314288  01 

0. 1083I428E  01 
0.108314288  01 

0. 10831428b  01 

U. 886999998  00 
0.886999998  00 

0.886999998  00 
0.8 8699999C  00 

0 . 88699999E  00 
0.886999991  00 

0. 88699999b  00 

0. 8471 1428L  00 
0.847114288  00 

0.847U4288'  00 
0.847114288  00 

0.B4711428E  00 
0.847114281  00 

0.847U428E  00 

-0.249000008-00 

-0.249000008-00 

-0.249000001-00 
-0. 249000001' -00 

-0.249000008-00 

-0.249000008-00 

-0 . 2  4900000  b -00 

0.293999991-01 

0.293999998-01 

0. 29399999 L-Ul 
0.293999991-01 

0. 29399999  8-rOl 
0.293999998-01 

0.293999996-01 

-0.303928578-00 

-0.303923578-00 

-0.303928578-00 

-0.303928571-00 

-0.303928571-00 

-0.303928571-00 

-0. 30392Hb7C-00 

0.30142857000 

-0.20833333001 

-0.604285711:  00 
-0.66260929E-16 

-0.447857141-00 
0. L4930555C-00 

0.  150000001-00 
0.49)432451-16 

0.47357143E-00 

-0.30310851L-16 

0.85714285E-01 
-0. 10416666C-00 

0. 

-0. 

-0. 74999999001 
-0.24671622E-L6 

-0.571428571-01 
0.21 L47105C-16 


-0. 

0.173611108-01 


0.18517999E  01 


0. 19880000E-00 


-0. 32449999E-00 


0.  1083142  8b  01 


0.836999991  00 


0.8471 1428E  00 


0. 2 4 90000 OE -00 


0.293999998-01 


0. 30392857000 
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checked  by  the  operator.2 

2' 

Because  NDEG  =  3,  after  storing  the  inverted  10  x. 

10  cubic  matrix,  the  program  proceeds  to  804  and  initializes 
for  Phase  3.  Now  the  program  is  ready  to  solve  Equation 
(ix)  for  the  polynomial  coefficients;  the d egree  1,  2  and 
3  coefficients  are  obtained  by  successive  operations,  'it 
851  the  linear  controls  are  read,  and  multiplication  of 
parts  of  the  stored  Y  and  B  matrices  develops  linear  coeffi¬ 
cients  in  the  matrix  C(I,J,K).  These  are  listed  in  floating 
point  notation  with  the  three  coefficients  (bo,  bp,  and  bg) 
in  columns  (Table  9);  successive  columns  relate  to  variables 
X^through  Xg(with  Xg  through  Xg  on  a  lower  line). 

Since  K  !S  1  the  machine,  after  listing  the  linear 
coefficients,  reads  the  degree  2  controls  at  852  and  returns 
to  699  to  compute  the  degree  2  coefficients  which  are  stored 
in  the  same  C(I,J,K)  matrix.  The  program  follows  the  same 
looping  so  that  the  six  coefficients  are  listed  in  columns 
(Table  10).  K  is  now  2,  so  the  program  goes  to  12  to 
test  whether  card  output  is  required. 

If  coefficient  cards  are  required  for  Part  II  of  the 
program  PCHOO  will  have  been  s  et  at  1  on  the  second  master 
card.  If  the  number  (NX)  of  dependent  variables  (Xn)  Is 
four  or  less  the  machine  proceeds  to  7  and  puts  Instructions 
on  Tape  7  for  the  punching  of  six  cards,  i.e.,  one  for  each 

coefficient  cq,  cp, . eg.  The  format  of  these  cards  insures 

that  j 

Columns  1  through  6  -  Project  number  of  program  (NPR0J1)„ 

Column  7  -  Degree  of  surface  where  2  indicates 

degree  2,  and  3  degree  3 

Column  8  -  1,  which  designates  that  coefficients 

of  group  Xp  through  X^  are  involved. 

Columns  9  and  10  -  Numbers  cards  of  each  degree  in  cor¬ 

rect  order,  successive  cards  being 
1  o  ^ 

2 

As  an  additional  check,  the  program  could  be  modified 

either  to 

(a)  invert  the  inverse  matrix,  A(I,J);  this  should 
recompute  the  original  matrix  P(I,J),  or 

(b)  multiply  the  original  matrix,  p(I,J),  by  the  com¬ 
puted  Inverse,  A(I,J);  the  product  should  be  a 
close  approach  to  the  unit  matrix. 


i 
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Columns  11  through  25  -  Coefficients  for  variable  Xi 

Columns  26  through  40  -  Coefficients  for  variable  Xg 

Columns  41  through  55  -  Coefficients  for  variable  X3 

Columns  56  through  70  -  Coefficients  for  variable  X4 

If  there  are  five  to  eight  variables  NX  is  greater 
than  4  and  the  program  goes  to  8,  instead  of  7,  so  that  two 
sets  of  coefficient  cards  are  punched.  The  first  set  con¬ 
tains  variables  X]_  through  X4  identified  by  1  in  column  8, 
and  the  second  contains  X5  through  Xq  and  is  identified  by  2 
in  column  8.  The  format  of  these  cards  is  otherwise  the  same, 
and  in  both  cases  the  coefficients  are  given  in  floating 
point . 


The  degree  2  controls  made  K  =  2,  so  the  program  pro¬ 
ceeds  to  853  and  reads  the  degree  3  controls  (including 
K  -  3).  The  degree  3  coefficients  are  computed  and  added  to 

C(I,J,K);  these  coefficients  (d^,  dg,  dg,  d4,  dg, . d10)  are 

listed  in  columns  (Table  11).  Cards  are  punched  when  PCHCO 
is  1,  with  the  format  detailed  above  for  the  degree  2. 

K  is  now  3  so  the  program  moves  to  854  where  FPHAS 
is  changed  to  2  and  phases  4  and  5  are  initiated.  Following 
initialization,  the  first  data  card  is  read  (at  25)  for  the 
second  time  (B  on  Flow  Sheet).  In  practice  it  is  the  first 
card  of  the  second  (duplicate)  data  deck  which  is  read.  Pro¬ 
ceeding  to  190,  the  data  cards  are  counted,  and  then  (re¬ 
turning  to  37)  powers  and  products  of  U  and  V  are  compiled 
for  this  first  data  card. 

The  linear  controls  required  for  computing  the  values 
of  Xn  on  the  trend  surface  at  each  U,V-grid  point  are  read  at 
200.  Each  dependent  observed  variable  (Xn)  is  stored  in 
SUMX(L)  and  as  subsequent  data  cards  are  read  the  values  of 
Xn  are  added  to  SUMX(L)  to  derive  the  XXi,  £X2, . ,  £Xs« 

Also,,  each  dependent  observed  variable  is  squared  and  stored 
in  SUMXQ(L)  to  develop  the  Ex^. 

At  710  a  brief  caption  stating  the  d egree  surface  and 
the  card  number  is  listed  in  preparation  for  tabulation  of 
the  main  answers. 

At  717  the  linear  computed  values  of  Xn  at  the  U,V- 
point  recorded,  on  the  first  data  card  are  calculated  by 
solution  of  Equation  (ii);  namely. 
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CUEFFS  DEGREE  3 

1 

O.10052714U  01 
0.L00527L4L  01 

0.10052714b  01 
0* 100527 14E  Oi 

2 

-U.3027S238E-00 

-0.30275238E-00 

-0*  3 02 752 38b -00 
-0. 302752386-00 

3 

0 . 241950  7Ih  01 

0. 24195071b  01 

0.24195071b  0.1 
0.241950716  01 

4 

0*  65300000b  00 
0.65300000E  00 

0. 65300000E  00 
0. 6  5300000E  00 

5 

0. 16547143fc-00 

0. 16547143E-00 

0.16547143000 

0.16347143000 

6 

-0.13942143b  01 
-0.13942143E  01 

-0.13942L43L  01 
-0. 13942143b  01 

7 

-0. 16733333E-00 
-0. 16733333t-00 

-0. 16733333000 
-0. 16733333000 

8 

-0.74499999C:-01 

-0.74499999E-01 

-0.  7449999900 1 
-0.74499999001 

9 

0.21S57142E-01 

0 . 2 1857  1.42t-0 1 

0.21857142001 
0. 2  1857 14200 1 

10 

0. 17625000E-00 

0. 1 7625000E-00 

0. 17625000E-00 
0. L  7625000t-00 

TABLE  11 


0. 100527 14E  01 

0. 100527 14E  01 

0.10052714b  01 

0.10052714b  01 

-0.30275238000 
-0. 30275238c- 00 

-0.302752386-00 

-0.302752386-00 

0. 24195071b  01 
0.241950716  01 

0.241950 71E  01 

0.2419  5071b  01. 

0.65300000c  00 
0.65300000b  00 

0.653000006  00 

0 • 65300000L  00 

0. 16547L43t-00 

0. 165471436-00 

0. 16547 1436-00 

0. 16547143t-00 

-0.13942143b  01 
“0.139421436  01 

-0.13942L43C  01 

-0.139421436  01 

-0.167333336-00 

-0.167333336-00 

-0. 167333336-00 

-0. L6733333L-00 

-0. 7  4499999C-0 1 
-0.744999996-01 

-0.744999996-01 

-0. 74499999E-01 

0.218571426-01 

0.218571426-01 

0.218571426-01 

0.218571426-01 

0.176250006-00 

0. 176250006-00 

0.176250006-00 

0.176250006-00 

TABLE  12 


CONTROL 


PT 

DEGREE,  l  CARD 

1 

X  OBS 

0001 

XI 

1.00000 

0001 

X2 

1.00000 

0001 

X  3 

l. 00000 

0001 

X4  * 

1.00000 

0001 

X5 

1.00000 

0001 

X6 

1.00000 

0001 

X7 

1.00000 

0001 

X8 

l. 00000 

DEGREE 

:  2  CARD 

1 

0001 

XI 

1.00000 

0001 

X2 

1.00000 

0001 

X3 

1.00000 

0001 

X4 

l. 00000 

0001 

X5 

1.00000 

0001 

X  6 

L. 00000 

0001 

X7 

1.00000 

0001 

X8 

1.00000 

DEGREE 

3  CARD 

1 

0001 

XI 

1.00000 

000.1 

X2 

l.  00000 

0001 

X3 

1.00000 

0001 

X4 

1.00000 

0001 

X5 

1.00000 

0001 

X6 

1.00000 

0001 

X7 

l. 00000 

0001 

X8 

1.00000 

DEGREE 

1  CARD 

2 

0002 

XI 

1.50000 

0002 

X2 

1.50000 

0002 

X3 

1.50000 

0002 

X4 

1.50000 

0002 

X5 

1.50000 

0002 

X6 

1.50000 

0002 

X  7 

1.50000 

0002 

X8 

1.50000 

DEGREE 

2  CARD 

2 

0002 

XI 

1.50000 

0002 

X2 

L . 50000 

0002 

X3 

1.50000 

0002 

X4 

1.50000 

0002 

X  5 

1.50000 

0002 

X6 

1*  500 00' 

0002 

X  7 

1.50000 

0002 

X8 

1.50000 

DEGREE 

3  CARD 

2 

0002 

XI 

1.50000 

0002 

X2 

1.50000 

0002 

X3 

1.50000' 

0002 

X4 

1.50000 

0002 

X5 

l .  50000 

0002 

X6 

L. 50000 

0002 

X7 

1 .  *>0000 

0002 

X'8 

L  .50000 

COMP 

DEVI  AT  ION 

1.85180 

-0. 85180 

1.85180 

-0. 85  L  80 

1.85180 

-0.8618C 

1.85180 

-0.85180 

1.85180 

-0.85180 

1.85180 

-0.85180 

1.85180 

-0.85180 

1.85180 

-0.8518U 

1.08314 

-0.08314 

1.08314 

-0.08314 

1.08314 

-0.08314 

1.08314 

-0.08314 

1.08314 

-0.08314 

1.08314 

-0.08314 

1.08314 

-0.08314 

1.08314 

-0.08314 

1.00527 

-0.00527 

1.00527 

-0.00527 

1,00527 

-0.00527 

1.00527 

-0.00527 

1.00527 

-0.00527 

1.00527 

-0.00527 

1,.  00527 

-0.00527 

1.00527 

-0.00527 

2.05060 

-0.55060 

2.05060 

-0.55060 

2.05060 

-0.55060 

2.05060 

-0.55060 

2.05060 

-0.55060 

2.05060 

-0.55060 

2.05060 

-0.55060 

2.05060 

-0.55060. 

1. 72114 

-0.22114 

1.72114 

-0.22114 

1.72114 

-0.22114 

1.72114 

-0.22114 

1.72114 

-0.22114 

1.72114 

-0.22114 

1.72114 

-0.22114 

1.72114 

-0.22114 

l.  LR819 

0.31181 

1  .  18819 

0.31181 

1.18819 

0.31181 

1.18919 

0 . 3 1 1  8  1 

l.  13819 

0. 31181 

L .  18819 

0 . 3  1  l  8  l 

1. 1 8819 

0  .  U  L  8  l 

l.  18319 

0 .  U  l  ft  ) 
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X  -  b0  +  bxU  +  b2V, 


with  the  aid  of  the  linear  coefficients  stored  in  C(I,J,K) 
and  the  powers  of  U  and  V  compiled  at  37  in  T(I).  The 
computed  values  are  stored  in  XCP(L,K).  These  computed 
values,  XCP(L,K),  are  subtracted  from  the  observed  values 
of  Xn  on  the  first  data  card.  The  differences  are  the 
deviations  at  this  U,V-point,  and  are  stored  in  RESD(L,K)» 

As  successive  data,  cards  are  read,  the  computed  values 
are  calculated  and  added  to  SXCP(L.,K)  to  provide  the  sum  of 
the  computed  Xn.  Similarly,  both  the  computed  values  and 
the  deviations  are  squared,  and  successive  squares  are  added 
to,  and  stored  in,  SXCPQ( L, K)  and  SRESQ(L,K),  respectively. 

The  principal  answer's  relating  to  the  first  datum- 
point  are  now  listed  in  fixed  point.  The  observed  value 
(X  OBS)  of  X-j_,  the  computed  Value  (X  COMP),  and  the  d  eviation 

are  listed  across  the  page,  and  the  values  for  X2,  X3, . 

follow  on  s uccessive . lines. 

The  linear  controls  at  200  defined  K  3  1  so  the  pro¬ 
gram  proceeds  to  862  where  the  degree  2  controls  are  read 
(including  K  a  2).  The  degree  2  equations ( iii)  are  now 
solved  with  values  from  the  first  data  card,  and  the  steps 
enumerated  above  are  repeated  to  develop  degree  2  values 
which  are  stored  in  the  same  matrices  with  the  linear  values., 
On  re-reaching  715  the  observed,  computed,  and  deviation 
values  of  Xn  are  listed  for  e ach  dependent  variable. 

The  program  loops  to  863  and  reads  the  degree  3  con¬ 
trols  (including  K  -  5)  and  returns  to  717  again  to  solve 
the  degree  3  polynomial  equation  for  values  of  Xn.  A  similar 
set  of  values  as  for  the  degree  1  and  2  surfaces  is  com¬ 
puted,  and  the  degree  3  answers  for  the  first  U,V-point 
are  listed  in  the  form  shown  in  Table  12. 

Since  K  =  3,  the  program  loops  back  to  25  (B  on  the 
Plow  Sheet)  to  read  the  second  data  card.  The  degree  1,  2 
and  3  answers  for  the  second  U,V-point  are  prepared  and  listed 
(Table  12),  and  the  new  values  are  added  to  the  matrices  being 
used  for  summation,  i.e.,  SUMX(E),  SUMXQ(L),  SXCP(L,K),  etc. 

It  will  be  noticed  that  the  sample  number  is  listed 
at  the  extreme  left  of  the  answer  tabulation. 

The  program  keeps  looping  back  to  25  until  all  the 
data  cards  have  been  read.  When  the  second  nines  card  is 
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TABLE  13 


SUMS  OF  SQUARES 

X  UBS 

X  COMP 

UEVIAJ ION 

DEGREE  L 

XL 

16,44038 

5.20005 

11.24033 

X2 

16.44038 

5.20005 

11.24033 

X'3 

16,44038 

5.20005 

11.24033 

X4 

16.44038 

5.20005 

11.24033 

X5 

16.44038 

5.20005 

1 1.24033 

X6 

16.44038 

5.20005 

11.24033 

X7 

16.44038 

5.20005 

11.24033 

X8 

16.44038 

5.20005 

11.24033 

DEGREE  2 

XL 

16.44038 

11-65615 

4.78 A 23 

X2 

16.44038 

11.65615 

4.78423 

X  3 

16.44038 

11.65615 

4.78423 

X4 

16.44038 

11.65615 

4.78423 

X5 

16.44038 

11.65615 

4.78423 

X6 

16.44038 

11.65615 

4.78423 

X7 

16.44038 

11.65615 

4.78423 

X8 

16.44030 

1  1.65615 

4.78423 

DEGREE  3 

XI 

16.44038 

13.95289 

2.48749 

X2 

16.44038 

1 3.95289 

2.48749 

X  3 

16.44038. 

13.95289 

2.48749 

X4 

16.44038 

13.95289 

2.48749 

X5 

16.44038 

13.95289 

2.48749 

X'6 

16.44038 

13.95289 

2.48749 

X7 

16.44038 

13.95289 

2.48749 

X8 

16.44038 

13.95289 

2.48749 

TABLE  14. 


•ERCONT  OF  SUM  OF 

SUJARES  ACCOUNTED 
L iNtAK 

FOR  BY  EACH  SURFACE 
QUAORAI IC 

CUBIC 

XI 

31.62972 

70.89951 

84.86964 

X2 

31.62972 

70.89951 

84.86964 

X  3 

31.62972 

70.89951 

84.86964 

X4 

31.6297 2 

70.89951 

84.86964 

X5 

31.62972 

70.89951 

84.86964 

X6 

31.02972 

70.89951 

84.86964 

X7 

31.62972 

70.89951 

84.86964 

X8 

3U62972 

70.89951 

84.86964 
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reached  at  the  end  of  the  deck,  the  nines  test  (180)  causes 
the  program  to  branch  to  255,  Titles  for  the  sums  of 
squares  summaries  are  now  listed  and  the  required  values 
are  computed  in  accordance  with  Equations  (xii)  through  (xv). 
These  values  are  stored  under  the  following  names: 


SSX(L) 
SSXCP(L, K) 

ssres(l,k) 

SSXCR(L,K) 


total  sum  of  squares  of  the  observed  data 
sum  of  squares  of  the  computed  values 
sum  of  squares  of  the  deviations 
percentage  of  total  variability  accounted 
for  by  each  trend  component. 


These  values  are  computed  for  each  dependent  variable,  and 
(except  in  the  case  of  SSX(L)  )  for  each  degree. 


On  reaching  718  the  first  .three  of  these  sets  of 
values  are  listed  across  the  page;  values  for  Xp,  X2,  Xg, 
are  written  on  successive  lines  (Table  15),  Degree  I 
(linear)  values  are  given  first  for  all  variables,  then 
degrees  2  and  3  (linear  plus  quadratic,  and  linear  plus 
quadratic  plus  cubic,  respectively). 


The  percentage  sums  of  squares  accounted  for  by  each 
surface  are  then  tabulated  in  fixed  point  notation  for  each 
degree  and  all  dependent  variables  (Table  14). 

7/hen  this  listing  is  completed,  the  program  loops 
back  to  3  and  is  ready  to  begin  the  next  completely  new 
problem.  There  is  no  limit  to  the  number  of  successive 
problems  which  can  be  run. 


PREPARATION  TO  RUN-  THE  PROGRAM 
(PART  II) 


The  data  deck: 

The  data  deck  is  assembled  in  the  following  manner: 

1.  vDATA 

2.  Pour  TITLE  CARDS 

3.  One  MASTER  SARD 

4.  Deck  of  data  cards 
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If  more  than  one  problem  is  to  be  executed  with  this 
program  in  the  same  run,  successive  decks  may  be  assembled 
according  to  2  through  4  above,  and  placed  behind  4  of  the 
first  deck..  Any  number  of  successive  problems  can  be 
executed  in  this  manner. 


Title  cards: 

The  format  is  (12A6)  so  that  numeric  and  alphabetic 
characters  can  be  utilized;  typical  cards  might  have  the 
following  form: 

1  WHITT  El?  PROJECT  040061 

0  MODAL  DATA  FOR  IACOHNE  PLUTON,  QUEBEC 

0  XI  -  QTZ;  X2  =  COLOR  INDEX;  X3  =  FELD.  RATIO 

0  LQC -TREND  SURFACES,  UV -COMPUTED  C-RID 


Master  card: 

This  is  a  complex  card  which  has  the  format  (A6,2F10,3, 

I3,2F10.3, 13, II)  and  it  must  provide  the  following  information: 

Columns  1  through  6  -  Project  control  number  (numeric  or 

alphabetic  characters)  identified  as 
NPROJ. 

Columns  7  through  16  -  The  maximum  U  coordinate  (UB  on 

Figure  9)  which  must  be  expressed  to 
three  decimal  places;  identified  as 
UB  in  the  program. 

Columns  17  through  26  -  The  minimum  U  coordinate  (UA  on 

Figure  9)  which  must  be  expressed  to 
three  decimal  places;  identified  as 
UA  in  the  program. 

Columns  27  through  29  -  The  inclusive  number  of  rows  (identi¬ 
fied  as  ITU)  of  computed  values  re¬ 
quired.  The  maximum  number  possible 
is  35 . 

Columns  30  through  39  -  The  maximum  V  coordinate  (VD  in 

program)  expressed  to  three  decimal 
places . 
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Figure  10:  Flow  chart  for  Part,  II  of  the  program 
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Columns  40  through  49  -  The  minimum  V  coordinate  (VC  in 

program)  expressed  to  three  decimal 
places . 

Columns  50  through  52  -  The  inclusive  number  of  columns 

(identified  as  NV)  of  the  computed 
values  required.  Maximum  number 
possible  is  35. 

Column  53  -  The  number  of  dependent  variables 

(Xn)  which  can  be  from1  one  to  four 
only. 

As  an  example,  to  obtain  computed  values  for  f|our 
dependent  variables  at  the  points  marked  with  a  X  in  ^igure 
9  the  master  card  must  have  the  form: 

040061000002 .000000000. 500004000002. 5 00000000. 5 000054  * 

Data  cards: 

These  comprise  the  entire  card  output  from  Part  I  of 
the  program.  When  from  one  to  four  dependent  variables  (Xn) 
were  involved  in  Part  I,  the  output  cards  comprise  the 
new  data  deck  without  modification.  This  deck  contains  six 
degree  2,  followed  by  ten  degree  3,  coefficient  cards.  If 
five  to  eight  dependent  variables  were  involved  in  Part  I, 

the  first,  third,  fifth, . .  and  thirty-first  cards,  which 

are  identified  by  1  in  column  8,  relate  to  variables  X]_ 

through  X4  while  the  second,  fourth,  sixth, . .  and 

thirty-second  cards,  which  are  identified  by  2  in  column  8, 
relate  to  variables  X5  through  Xq  .  In  Part  II  it  is  neces¬ 
sary  to  run  these  two  groups  as  successive  separate  problems 
(the  master  and  title  cards  can  of  course  be  identical  for 
both  groups),  because  it  will  only  accommodate  four  dependent 
variables  at  a  time. 

No  nines  cards  are  required. 


OPERATION  OF  THE  PROGRAM 
(  PART  II ) 

The  sequence  of  operations  can  be  read  from  the 
flow  sheet  (Figure  10)  and  is  described,  below.  In  this 
illustration,  the  card  output  from  the  Illustration  in  Part 
I  is  assumed. 
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TABLE  15 

r 

WHITTEN  PROJECT  040108 

UV-COMPUTED  VALUES  ^ 

XL  THROUGH  X4  IDENTICAL  SYNTHETIC 
PROBLEM  TO  ILLUSTRATE  PROGRAM 


TABLE  16 


V-COOROt NATES 


0. 

1.000 

2.000 

3.000 

4.000 


U-COORD I  NAVE  S 


0. 

1.000 

2.000 

3.000 


TABLE.  17 


UV-COMPUTEl)’  VALUES  LISTED  BY  ROWS  ACROSS  MAP 
QUADRATIC  ANSWERS  FIRST  FOLLOWED  BY  CUBIC 

VARIABLE  X-L 


QUADRATIC  COMPUTED 

VALUES 

u= 

0. 

1.083 

1.626 

1.562 

0.889 

-0.391 

u= 

UOOO 

1.721 

2. 294 

2.258 

1.615 

0.364 

u= 

2.000 

1.861 

2.46  3 

2.457 

1.844 

0.622 

u= 

3.000 

1.503 

2.  135 

2f  158 

1.574 

0.382 

CUBIC 

COMPUTED  VALUES 

U- 

0. 

1.005 

2.207 

1.677 

0.475 

-0.344 

U=' 

1.000 

1.  188 

2.503 

2.130 

1.  127 

0.552 

u= 

2.000 

1.673 

2.95  L 

2.586 

1.635 

1.155 

u= 

3.000 

1.456 

2.549 

2.042 

0.993 

0.459 

VARIABLE  X- 

2 

QUADRATIC  COMPUTED 

VALUES 

u= 

0. 

1.083 

1.626 

1.562 

0.889 

-0.391 
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The  four  title  cards  are  read  first  and  this  informa¬ 
tion  is  immediately  printed  out  (Table  15). 

The  master  card  is  read  and  the  information  stored. 
After  initialization,  six  degree  2  coefficient  cards  (B  in 
the  flow  sheet),  and  the  ten  degree  3  coefficient  cards  (C 
in.  the  flow  sheet),,  are  read  and  stored  in  matrix  COE(I,J). 

The  U  and  V-coordinates  for  the  required  map-area  are 
determined  next.  The  interval  between  successive  rows 

UU  =  (UB  -  UA)  /  (UN  -  1) 
and  successive  columns 

VV  =  (VD  -  VC)  /  (VN  -  1) 

are  computed  from  information  on  the  master  card.  These  in¬ 
tervals  are  successively  added  to  UA  and  VC  and  stored  in 
UUM(I)  and  VVM(I),  respectively.  At  103  a  list  of  the  V- 
coordinates,  VVM(I),  and  of  the  U-coordinates,  UUM(I),  is 
printed  out  (Table  16). 

The  captions  for  the  computed  grid  of  values  is 
printed  at  107.  Then  M  (set  as  1)  defines  which  variable 
is  involved  in  the  following  calculations.  Using  the 
coefficients  stored  in  C0E(I,J),  and  the  U  and  V~coordinates 
stored  in  UUM(I)  and  VVM(I),  respectively,  the  degree  2  and 
degree  3  polynomials  are  solved  with  respect  to  Xq  for  each 
U,V-grid  point.  The  values  are  stored  in  C0MPX(I, J,K) . 

This  operation  is  effected  by  substituting  appropriate  values 
in  Equation  (i).  The  degree  2  values  are  then  listed  in 
fixed  point  to  three  decimal  places  by  rows  across  the  page, 
and  the  U-value  of  the  row  is  identified  in  the  left  margin. 
The  degree  3  values  for  Xp  are  then  listed  in  a  similar  manner 
(Table  17). 

Following  this  listing,  the  program  compares  the 
number  of  variables  (NVAR)  stated  on  the  master  card  with 
M.  If  there  is  still  another  variable,  the  program  incre¬ 
ments  M  to  ( M  + 1 ) ,  loops  back  to  19,  and  initializes  the 
C0MPX(I,J,K)  matrix.  The  coefficients  for  Xg  are  then  used 
to  compute  an  array  of.  values  which  are  stored  in  C0MPX(I,J,K) 
and  listed.  This  process  is  repeated  -until  values  for  all 
the  Xn  involved  have  been  computed  and  listed.  The  program 
then  loops  back  to  3333  to  begin  the  next  complete  problem. 
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APPENDIX  I 


NORTHWESTERN  UNIVERSITY  GEOLOGY  DEPARTMENT 

C  SURFACE-FITTING  PROGRAM  FOR  IRREGULARLY-SPACED  MAP  DATA 

C  LINEAR,  QUADRATIC,  CUBIC  -  DOUBLE  PRECISION  -  EIGHT  VARIABLES 

C  1 957  KRUMBEIN  AND  HARRIS  MACHINE  LANGUAGE  FOR  BASIC  IBM  650  4  VAR 

C  1 960  KRUMBEIN  AND  FAULKNER  REWRITTEN  IN  SOAP  FOR  IBM  650  -  4  VARS', 

C  1961  AXELROD  AND  BENSON  FOR  KRUMBEIN  REWRITTEN  IN  FORTRAN  FOR 

C  IBM  709  -  FOUR  VARIABLES  AND  SINGLE  PRECISION 

C  1962  WHITTEN  EXTENDED  TO  DOUBLE  PRECISION  EIGHT  VARIABLES  AND 

C  CARO  OUTPUT  FOR  COEFFICIENTS  FOR  COMPUTING  GRID  OF  UV-VALUES 

DIMENSION  ZERO  {  2  )  ,  0NEI.2) 

D  DIMENSION  U ( L )  ,  VII) 

DIMENSION  T I  T  l  (  9  )  ,  TIT2(9),  T IT3 ( 9 ) , T  l  T4 ( 9 1 , VF T{ 1 0 ) 

D  DIMENSION  X  (  8  )  ,  TUO),  B  { 1 0  *  8  ) 

0  DIMENSION  Y(10,10,3),  C(10,8,3) 

DIMENSION  XCP<8,3),  RESD(8,3) 

DIMENSION  SUMXI8),  SUMXQI8),  SSXI8) 

DIMENSION  SXC P { 8 , 3 )  ,  SXCPC ( 8 , 3  )  ,  SSXCP(8,3) 

DIMENSION  SRE SQ ( 8,3),  SSRES ( 8 , 3 ) , SSXCR ( 8, 3 ) 

D  DIMENSION  A ( 1 0 , 10 ) , P ( 10 , 1 0 ) 

CARO  1=  LH l 
CARD2  =  IH2 
ONE  I i  )  =  1.0 
ONE ( 2  )  =  0.0 

3  READ  INPUT  TAPE  5 , 72 , T I T L , T I T2 , T I T 3 , T I T4 
WRITE  UUTPUT  TAPE  6 , 72 , T IT l , T I T2 , T I T3 , T I T4 

72  FORMAT  (9A6) 

C  *  INITIALIZE  PHASE  ONE 

ZEROU)  =  0.0 
ZERO  (2)  =  0.0 

5  NPHA S=  l 

DO  20'  1  =  1,  10 
DO  10  J=l,10 
0  10  PI  I,J)  =  ZERO 

DO  20  J  =  1,8 
D  20  8 ( I , J  )  =  ZERO 

DO  904  I  =  1,10 

D  904  T(I)  =  ZERO 
D  U  =  ZERO 

D  V  =  ZERO 

0  TCI)— 1.0 

C  READ  MASTER  CARDS 

READ  INPUT  TAPE  5 , 73 , NPRO J l , NPR0J2 , NX , PCHCO 

73  FORMAT  ( Afa , A4 ,  1 1 , F  L . 0 ) 

READ  INPUT  TAPE  5,li,VFT 

11  FORMAT  ( 1 OX , 1 0  A  6 ) 

WRITE  OUTPUT  TAPE  6,,  9  ,  (  I  ,  I  =  1 ,  NX  ) 

9  FORMAT ( /20X4HDATA/20X4H - //5X 7HCONTROL , 4X IHU, 10X 1HV , X , 8 { 9X IHX , I 1 

l  )  ) 

C  *  ADO  IN  DATA 

25  READ  INPUT  TAPE  5 , VF T , I D , U , V , ( X (  I  ) ,  I  =  1 , NX ) 

GO  TO  ( 27,  180)  , NPH AS 

27  IF  (U-99999.)  30,31,31 

30  CONTINUE 

WRITE  OUTPUT  TAPE  6,  80 1 3 , I D , U , V , ( X (  I  )  , I  =  1 , NX ) 

8013  FORMAT  <5X,  A4,10Fil.3) 

C  ASSEMBLE  MATRIX 


o  O  OOOOO  OO  O  o  CTOO  oooooocoo 
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37  T ( 2 ) =U  I 

T ( 3 ) =V  I 

T(4)=U*U 

r  <  5 ) =u* v 

T ( 6 ) = V* V 

T ( 7 ) =U*U*U 

T ( 8 ) =U*U* V 

U9)=u»v*v 

T ( 10)=V*V*V 

GO  ro  < 38,200) ,NPHAS 

38  DO  45  1  =  1 ,  10 
DO  40  J=l,10 

40  P{I,J)  =  T(I)*TU)+P(I,J> 

DO  45  J  =  1 » NX 

45  B(I,J)=X(J)*T(I)+8(I»J) 

GO  TO  25 

WRITE  UV  CUBIC  MATRIX 

31  WRITE  OUTPUT  TAPE  6,78, 

78  FORMAT  ( /// /10X8HUV  CUBIC/) 

DO  32  1=1,10 

32  WRITE  OUTPUT  TAPE  6 , 70 , I , ( P ( I , J ) , J= 1 , 10 ) 
70  FORMAT  ( / 3X I  2 , 5E20 . 8/ ( 5X , 5E20 . 8 ) ) 

WRITE  OUTPUT  TAPE  6,79, 

WRITE  OUT  VECTOR  XS 

79  FORMAT  ( / /  10X1 1HVEC TOR  X IS)/) 

DO  33  1=1,10 

33  WRITE  OUTPUT  TAPE  6 , 70 , I , ( 8 ( I , J ) , J= 1 , MX ) 
BEGINNING  PHASE  TWO 

00  501  1=1,10 
DO  501  J=  1 » 10 
00  501  K=  1 ,3 
501  Y ( I , J , K )  =  ZERO 

*  CONTROL  FOR  INVERSION  OF  DATA  MATRIX 
WRITE  OUTPUT  TAPE  6,74, 

801  MH=3 

NOEG  =1 
GO  TO  600 

802  MH=6 
NDEG  =2 
GO  TO  600 

803  MH= 1 0 
NDEG  =3 
GO  TO  600 

600  DO  514  1=1, MH 
DO  514  J= 1 » MH 
514  A ( I,J)  =  P(  I.J) 

MATRIX  INVERSION 
DO  414  K= 1 , MH 
DIV  =  A ( K , K ) 

A ( K , K  )  =  ONE 
DO  411  J=1 , MH 

411  A(K,J)=A(K,J)/OIV 
DO  414  I = 1 , MH 

IF! I-K)412,414,412 

412  DIV  =  All ,K) 

A! I ,K)  =  ZERO 
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00  413  J=1,MH 

D  413  A( I , J )=A( I , J)-0IV*A(K, J) 

414  CONTINUE 

WRITE  OUTPUT  TAPE  6,900,NDEG 

900  FURMAI  (/10X.26H INVER  TED  MATRIX  OF  DEGREE  ,11) 

DO  140  1=1, MH 

140  WRITE  OUTPUT  TAPE  6 , 70 , I , ( A (I , J ) , J= l , MH ) 

DO  601  I = 1 , MH 
DO  601  J=1 , MH 
D  601  Y( I, J.NDEO!  =A(  I , J ) 

GO  TO  (802,803,804)  ,NDEG 
C  INITIALIZE  PHASE  THREE 

804  DO  300  K= 1 , 3 
DO  300  J- 1 , NX 
DO  300  1=1,10 
D  300  C { I , J , K )  =  ZERO 
C  COMPUTE  CGEFF 

WRITE  OUTPUT  IAPE  6,74, 

74  FORMAT  ( 1H1 ) 

851  K=  1 
MH  =  3 

GO  TO  699 

852  K  =  2 
MH  =  6 

GO  TO  699 

853  l<  =  3 
MH=  1 0 

GO  TO  699 
699  DO  700  L= 1 , NX 
00  700  1=1, MH 
DO  700  J= 1 , MH 

D  700  C ( I , L  ,  K )  =  C ( I >  L , K )  +  Y(I,J,K)  *  B(J,L) 

WRITE  OUTPUT  (APE  6,901,K 

901  FORMAT  (/iOX,  14HC0EFF S  DEGREE  ,11) 

DO  703  1=1, MH 

703  WRITE  OUTPUT  IAPE  6 ,70 , 1 , t C (  I  , L , K )  ,  L=1,NX> 

GO  TO  (2, 12, 12), K 
12  IF  ( PCHCO )  1,2,1 

1  IF  (NX-4)  7,7,8 

7  DO  9001  I  = 1 , MH 

9001  WRITE  OUTPUT  TAPE  7 , 9000 , NPROJ 1 , K , CARD1 , I , ( C (  1 , L , K > , L=  1 , NX  ) 

9000  FORMAT  (A6, I  1 , A1 , 12,4015.8) 

GO  TO  2 

8  DO  9002  1=1, MH 

9002  WRITE  OUTPUT  TAPE  7 , 9000 , NPROJ 1 , K , CAR01 , l , ( C ( I , L , K ) , L= 1 , 4 ) , 

1  NPROJ l ,K,CARD2, 1 , (C( I , L »  K ) , L  =  5, NX ) 

2  GO  TO  (852,853,854)  ,  K 

854  NPHAS  =  2 

IF  (PCHCO)  4,6,4 
4  END  FILE  7 

C  BEGIN  COMPUTING  ANSWERS  -  PHASES  FOUR  AND  FIVE 

6  WRITE  OUTPUT  TAPE  6,75, 

75  FORMAT  (  8H ICON  I ROL / 3X2HPT  ,  15X5HX  0BS.14X6HX  COMP ,  13X9HDEV  I  A T  ION/ /  ) 
KCD  =  0 

DO  705  L= 1 » NX 
DO  704  K.=  1 , 3 
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SXCP (  L  »  K  )  =  0.0 
SXCPQ ( L  * K )  =0.0 

704  SRESQ(L.K)  =0.0 

SUMX(L)  =0.0 

705  SUMXQ(L)  =0.0 

C  *  READ  DATA 

GO  TO  25 

180  IF  (U  -  90999.)  190,255,255 
190  KCD=KCD+ 1 
GO  TO  37 

C  *  COMPUTE  X 

200  K= 1 
MH  =  3 

00  201  L= 1 ,NX 
SUMX(L)=SUMX(L)+X(L) 

201  SUMXQ(L)=SUMXQ(L)+X(L)*»2 
GO  TO  710 

862  K=2 
MH=6 

GO  TO  710 

863  K=3 
MH=  1 0 

GO  TO  710 

710  WRITE  OUTPUT  TAPE  6,902,K,KCD 

902  FORMAT  ( 8H  DEGREE  , 1 1 , 6H  CARD  ,13) 

DO  715  1=1, NX 

XCP ( L , K )  =0.0 
DO  717  1=1, MH 

717  XCP  ! L , K )  =C ( I »  L , K )  *  T ( I )  +  XCP(L.K) 

RESD  (L,K)  =X ( L ) -XCP ( L , K ) 

SXCP(L.K)  =  SXCPIL.K)  +  XCP ( L , K ) 

SXCPQ ( L , K )  =  S  XCPQ ( L , K )  +  XCP(L,K)**2 
SRESQ ( L , K )  =  SRESQ ( L ,K )  +  RESD(L,K)**2 

715  WRITE  OUTPUT  TAPE  6, 71, ID, L,  X{L),  XCP(L,K),  RESD(L.K) 

71  FORMAT  (2XA4.3H  X,I1,3F20.5) 

GO  TO  (862,863,864)  ,K 

864  GO  TO  25 

255  WRITE  OUTPUT  TAPE  6,76, 

76  FORMAT  ( IH1 , 15HSUMS  OF  SQUARES/20X5HX  OBS,  14X,6HX  COMP  1 3X9HDE V I  AT 
1  ION// ) 

CD=KCD 

DO  718  K=  1  >  3 

WRITE  OUTPUT  TAPE  6,903, K 

903  FORMAT  (/8H  DEGREE  ,11) 

DO  718  L= 1 , NX 

SSX ( L )  =SUMXQ(L )-( SUMX (L) *»2/CD) 

SSXC  P ( L , K )  =  SXC  PQ ( L , K )  -  { SXCP ( L , K ) **2/CD ) 

S  SRE  S ( L , K )  =  SRESQ ( L ,K ) 

S SXCR  ( L , K  )  =  ( SSXCP ( L , K ) /SSX ( L ) ) *100.0 

718  WRITE  OUTPUT  TAPE  6,77  , L , SSX (L ), SSXCP ( L , K ), SSRES ( L , K ) 

WRITE  OUTPUT  TAPE  6,9078, 

9078  FORMAT  (///2X»55H PERCENT  OF  SUM  OF  SQUARES  ACCOUNTED  FOR  BY  EACH 
9078  l  SURFA.CE/20X ,  6HL  INEAR,  1 3X,  9HQUADRA  T  IC  ,  10X  ,  5HCUB  IC ) 

DO  97.18  L  =  l ,  NX 

9718  WRITE  OUTPUT  TAPE  6 , 77 , L , ( SS XCR ( L , K ) , X  =  1,3) 

77  FORMAT  (6X3H  X,I1,3F20.5) 

GO  TO  3 

END( l ,0,0, 0,0, 0,0, 0,0, 1,0,0,0,0,01 
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APPENDIX  II 


NORTHW'E  STERN  UNIVERSITY  GEOLOGY  DEPARTMENT 

C  GRID  OF  COMPUTED  VALUES  FROM  COEFFICIENTS  OF  LQC  PROGRAM 

C  PREPARED  bY  WHITTEN  -  JANUARY  1962 

DIMENSION  UUMJ35),  VVM (35 ) , COMPX ( 35 , 3 5 , L 1 ) 

DIMENSION  T I T l t I 2 ) »  TIT2I12),  TIT3(I2),  TIT4I12),  COE(20,4) 

C  *  *RE AD  TITLE  CARDS 

3333  READ  INPUT  TAPE  5 , 101 , TI T1 , T I  I  2 , TI T3 , TI T4 

WRITE  OUTPUr  TAPE  6,101,TITl,  TIT2,  TIT3,  TIT4 
101  FORMAT  (12A6! 

C  **  READ  MASTER  CARD 

READ  INPUT  TAPE  5,1,  NPROJ , UB , UA , NU , VO , VC , NV , NV AR 

1  FORMA!  (A6.2FT0.3, I3.2F10.3, 13, I  1) 

C  **  READ  QUADRATIC  COEFFICIENT  CARDS 

DO  10  1-=  1,20 
DO  10  J  =  l , N V AR 

10  CGE  (  I  ,  J  )  =  0-0 
DO  11  I  =  1,6 

11  READ  INPUT  TAPE  5,2,  (COE  (I,J>,  J  =  l.NVAR) 

2  FORMA  I  (10X,  4E15.8) 

C  **  READ  CUBIC  COEFFICENT  CARDS 

DO  13  I  =  11,20 

13  READ  INPUT  TAPE  5,3,  (COE  (I,J>,  J  -  l.NVAR) 

3  FORMA  I  ( 10X,  4C15.8! 

C  **  COMPUTE  UV-MATRIX 

UN  =  NU 
VN  •=  NV 

UU  =  (UB  -  UA ) / ( UN  -  1.0) 

VV  =  (VO  -  VC ) / ( VN  -  1.0) 

DO  14  I  =  1,  NU 

14  UUM(I)  =  UA  +  ( IJU*FL0ATF (  I  -1)) 

DO  L  5  I  =  1,  NV 

15  VVM(I)  =  VC  +  ( VV*FLOATF( I-l ) ) 

WRITE  OUTPUT  TAPE  6,103, 

103  FORMAT  (//5X,  l 3HV-C00RD I NA TCS ) 

WRITE  OUTPUT  TAPE  6 , 104 , ( VVM ( I ) ,  I  =  1,NV) 

104  FORMAT  (//(5X.F7.3)  ) 

WRITE  OUTPUT  TAPE  6,105, 

105  FORMAT  (//5X, 13HU-C00R0I NATES) 

WRITE  OUTPUT  TAPE  6 , 106 , ( UUM ( 1) ,  I  =  1,NU> 

106  FORMAT  (//(5X.F7.3)  ) 

WRITE  OUTPUT  TAPE  6,107, 

107  FORMAT ( / /5X , 44HU V-COMPUTE D  VALUES  LISTEO  BY  ROWS  ACROSS  MAP//5X, 
10714 1HQUADRAT IC  ANSWERS  FIRST  FOLLOWED  BY  CUBIC) 

M  =  l 

19  DU  12  K  =  1,11,10 
DO  12  J  =  1 , NV 
DU  12  I  =  1,  NU 

12  COMPX  ( I , J , K )  =  0.0 
DO  16  K  =  1,11,10 
DO  16  J  = 1 , NV 

DO  16  I  =  1 , NU 

16  COMPX ( I , J , K )  =  COE  ( K , M )  +1C0E  ( K + l , M ) * UUM ( I ) ) + ( C OE  ( K+2, M ) *VVM( J ) 
16.1  )  +  (  COE  (K  +  3,M)*UUM(  I  )**2)  +  (.C0fc  (  K  *4,  M  )  *UUM  (  I  )  *VVM  (  J  )  )  +  (  COE  (K  +  5,M) 
162*  VVM  (  J  )  *  *2  )  +  (  COE  (  K  +  6  ,  M )  *IJUM  (  I  )  *  *3 )  +  (  (COG  (  K  +  7 ,  M  )  *UUM  (  I  )  **2  ) 
163*VVM( J ) ) + ( (COE  ( K+8 , M ) *UUM ( I ) ) * VVM ( J ) * *2 ) + (COE  ( K+9, M ) *VVM ( J ) **3) 

WRITE  OUTPUT  TAPE  6,110,M 


56 


1L0  FORMAT  ( ///5X,  1 IHVARI  ABLE  X-,U) 

DO  1000  K  =1,11.10 
IF  (l<  -  l)  111,111,112 

111  WRITE  OUTPUT  IAPE  6,113, 

113  -FORMA  I  I // 10X , 2 5H QUADRAT I C  COMPUTED  VALUES//) 

GO  TO  99 

112  WRITE  OUTPUT  I  APE  6,114, 

114  FORMAT  ( // 10X , 2 1HCUBI C  COMPUTED  VALUES/) 

99  DO  1000  I  =1 , NU 

1000  WRI  TE  OUTPUT  TAPE  6*  102 , UUM(  I  )  ,  (  COMPX.I  I ,  J , K  )  ,  J=  1  ,MV  ) 
102  FORMAT  ( //2X , 2HU= , F7. 3 , 9F 1 2 . 3/ ( 1 l X , OF l 2 . 3 ) ) 

IF  (M  -  NVAR )  17,3333,3333 
17  M  =  M  +  1 
GO  TO  19 

END (1,0, 0,0, 0,0, 0,0, 0,1, 0,0, 0,0,0) 
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